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ABSTRACT  . 


Spectral  analysis  provides  powerful  techniques  for  the 
approximate  linear  analysis  of  the  input-output  relations  of  a  system 
When  a  system  is  subjected  to  a  stochastic  ,  band-limifed  driving 
function  estimates  of  the  frequency  response  function  ,  the  coherence 
function  and  the  information  transmission  rate  may  be  obtained  from 
consistent  spectral  and  cross-  spectral  estimates  , 

One  problem  in  applying  spectral  analysis  to  a  primary  sensory 
neurone  is  that  the  response  of  the  neurone  is  a  spike  train  .  A 
method  of  obtaining  regular  ,  alias-free  samples  of  a  spike  train 
with  an  arbitrary  sampling  rate  has  been  developed  ,  Using  a  package 
for  consistent  spectral  and  cross-spectral  estimation  the  frequency 
response  functions  ,  coherence  functions  and  information  transmission 
rates  of  two  insect  mechanoreceptors  have  been  estimated  »  Both 
these  receptors  have  frequency  response  functions  indicating  a  high- 
pass  filter  characteristic  ,  and  a  coherence  function  which  is  never 
greater  than  0.5  .  The  low  coherence  function  indicates  that  the 
response  of  the  receptor  is  noisy  and  nonlinear  ,  and  so  the  frequency 
response  function  is  not  an  accurate  characterization  of  the  input- 
output  characteristics  of  the  receptors  . 
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CHAPTER  1  . 


INTRODUCTION  . 

Following  Pringle  and  Wilson’s  (1952  )  use  of  linear  systems 
theory  to  characterize  the  input-output  relations  of  an  insect 
mechanoreceptor  linear  systems  theory  has  been  used  to  investigate 
the  signal  transfer  characteristics  of  a  large  number  of  single 
neurones  .  It  has  been  applied  to  receptors  (  Chapman  and  Smith, 

1963  ;  Borsellino,  Poppele  and  Terzuolo,  1965  ;  Brown  and  Stein,  1966 
Pinter,  1966  ;  Houk  and  Simon,  1967  ;  Dodge,  Knight  and  Toyoda,  1968 
Matthews  and  Stein,  1969a;  Poppele  and  Bowman,  197o  ;  Knight,  Toyoda 
and  Dodge,  1970  ;  Pearson  and  Holden,  1970  ;  Holden  and  French,  1971a 
1971b  )  ,  higher  order  sensory  neurones  (  Spekreijse,  1969  ;  Maffei, 
Cervetto  and  Fiorentini,  1970  ;  Melville  Jones  and  Milsum,  1970  )  , 
neuro-muscular  effectors  (  Partridge,  1965  ;  McKean  et  al.  ,  1970  )  , 
and  synapses  (  Hagiwara  and  Morita,  1962  ;  Watanabe,  1968  )  ,  Some 
of  these  ,  and  other  examples  ,  are  reviewed  in  the  Neurosciences 
Research  Program  Bulletin  on  ’Neural  Coding’  (  Perkel  and  Bullock, 
1968  )  and  the  Brainerd  conference  on  ’Systems  Analysis  in  Neuro¬ 
physiology’  (  Terzuolo,  1969  )  .  Linear  systems  theory  provides  a 
powerful  body  of  techniques  for  investigating  the  input-output 
relations  of  linear  systems  ,  and  can  also  be  used  as  a  general 
approach  to  the  approximate  analysis  of  nonlinear  systems  , 

In  spite  of  its  widespread  application  in  the  study  of  single 
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units  in  the  nervous  system  there  are  problems  in  treating  a  neuronal 
element  as  a  system  with  a  clearly  defined  input  function  x(t)  , 
or  set  of  input  functions  ,  and  a  clearly  defined  output  y(t)  .  For 
a  primary  sensory  neurone  ,  the  input  is  a  continuous  time  varying 
function  which  may  also  vary  spatially  ,  and  the  output  may  be  a 
continuous  function  ,  as  in  photo-receptors  (  Arden,  1969  )  or  a  crab 
mechanoreceptor  (  Bush  and  Roberts,  1968  )  ,  but  is  generally  a  spike 
train  .  For  a  central  neurone  ,  the  input  may  be  taken  as  a  large 
number  of  spike  trains  ,  and  the  output  as  a  single  spike  train  .  For 
a  neuro-muscular  effector  ,  the  input  is  a  spike  train  ,  or  number 
of  spike  trains  ,  and  the  output  is  a  continuous  function  (  say 
tension  )  .  Thus  in  most  cases  one  or  more  of  the  signals  is  a  spike 
train  .  Although  a  recorded  spike  train  is  a  continuous  band-limited 
signal  ,  since  the  duration  of  an  action  potential  is  short  compared 
to  the  mean  inter-spike  interval  ,  and  action  potentials  are  generally 
all-or  none  ,  it  is  usual  to  treat  a  spike  train  as  a  series  of 
identical  events  .  An  exception  to  this  is  given  in  Gestri,  Maffei 
and  Petracchi  (  1967  )  and  Barbi  and  Petracchi  (  1971  )  . 

If  all  action  potentials  in  a  spike  train  are  taken  as  equivalent 
events  there  are  two  ways  in  which  a  spike  train  may  be  treated  : 

i) .  a  continuous  function  may  be  generated  from  the  spike  train 
(  Terzuolo,  1970  )  .  This  may  be  the  averaged  interval  between 
spikes  ,  the  rate  of  spikes  averaged  over  some  time  or  number  of 
spikes  ,  or  an  estimate  of  the  probability  density  of  spikes  relative 
to  some  time  marker  (  the  cycle  or  post-stimulus  time  histogram  ). 
Alternatively  ,  unaveraged  ,  discontinuous  functions  may  be  generated  . 
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ii) .  the  spike  train  can  be  considered  as  a  realization  of  a 
stochastic  point- process  (  Moore,  Perkel  and  Segundo,  1966  ;  Stein, 
1970  )  or  as  a  series  of  unit  impulses  or  Dirac  delta  functions 
(  French  and  Holden,  1971a  ;  Knox,  1969  ) 

When  a  suitable  form  of  describing  the  spike  train  has  been 
chosen  linear  systems  theory  can  be  used  to  derive  a  relationship 
between  the  input  signal  and  the  output  signal  .  Such  an  approach 
is  purely  descriptive  -  the  neuronal  element  is  treated  as  a  ’black 
box'  and  a  mathematical  model  is  formulated  which  attempts  to 
describe  the  input-output  relations  of  the  black  box  .  However  ,  such 
a  quantitative  description  is  of  value  in  its  own  right  ,  and  may 
be  used  as  a  basis  for  either  dissecting  the  black  box  into  possible 
components  (  see  ,  for  example  ,  Terzuolo  et  al.  ,  1968  )  or  for 
building  models  of  simple  neuronal  networks  . 

A  variety  of  methods  ,  both  in  the  time  domain  and  in  the 
frequency  domain  ,  are  available  for  obtaining  a  quantitative 
description  of  the  input-output  relations  of  a  system  .  These  methods 
are  reviewed  in  Chapter  2  and  their  implementation  discussed  in 
Chapter  3  .  In  later  chapters  some  of  these  methods  are  applied 
to  two  sensory  receptors  ,  the  tactile  spine  of  the  cockroach  and  the 
chordotonal  organ  of  the  locust  . 


' 
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CHAPTER  2  . 


LINEAR  SYSTEMS  THEORY  . 

2.1.  Definition  of  a  linear  system  . 

A  linear  transformation  from  a  vector  space  to  a  vector  space 
is  a  function  A  which  provides  a  one-to-one  association  with 
each  vector  z  in  Z^  with  a  vector  A(z)  in  Z^  such  that: 


A(  +  z2  )  =  A(z1)  4-  Mz2) 


and 


2.1 


A(az) 


aA(z) 


for  all  vectors  z  in  Z^  and  all  scalars  a  .  If  the  vector  space 
Z^  contains  all  possible  input  functions  z  to  a  system  and  the 
vector  space  Z^  contains  all  possible  responses  of  the  system  then 
the  system  is  linear  if  and  only  if  equation  2.1  is  valid  for  all 
input  functions  z  and  all  responses  A(z)  .  This  is  a  statement  of 
the  principle  of  superposition  .  If  the  principle  of  superposition 
holds  for  only  some  of  the  possible  inputs  ,  and  not  for  other  inputs  , 
then  the  system  is  not  linear  . 


2_.  2 . _ Application  to  linear  ,  noise-free  systems  . 

For  a  linear  noise-free  system  with  a  single  input  or  driving 
function  x(t)  and  a  single  output  or  response  y(t)  the  output  is 
related  to  the  input  by  the  convolution  integral  : 


y(t) 


2.2 


x(t)  g(t,x)  dx 


and  if  the  system  is  also  time  invariant  a  function  g(t)  exists 
such  that 

y (t)  =  /q  x(t)  g(t  -t)  dx  2.3 


where  g(t)  is  the  weighting  function  of  the  system  . 

Treatment  of  equation  2.3  is  greatly  simplified  if  it  is  examined 
in  the  s  domain  ,  where  s  is  a  complex  variable  .  This  may  be 
achieved  by  use  of  the  Laplace  Tranform  (  Smith,  1966  )  .  The  Laplace 
Transform  is  a  linear  transform  and  is  defined  by  : 


L.T.  {y (t) } 


Y(s) 


00 


-st 

y  e 


dt 


2.4 


and  its  inverse  is  given  by 

L.T.'1  {Y (s) }  =  ^  /C+L  Y(s)  etSds  2.5 

2tt  c— j 00 

where  c  is  a  constant  which  must  be  defined  to  ensure  convergence 
in  this  contour  integration  .  Thus  y(t)  and  Y(s)  form  a  Laplace 
Transform  pair 

y  ( t )  Y(s) 


In  practice  it  is  generally  not  necessary  to  evaluate  the  forward 
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and  inverse  transforms  as  they  are  available  in  tables  for  most 
commonly  encountered  functions  . 

The  most  useful  properties  of  the  Laplace  Transform  are: 

a)  Linearity 

k  y(t)  ++  k  Y(s) 

y  (t)  +  f(t)  ■<—>■  Y(s)  +  F  (s) 

b)  Time  translation 

y(  t  -t  )  e”XSY(s) 

c)  Real  convolution 

/q  y(  t  -t  )  f  (t)  Y(s)F (s) 

Application  of  the  convolution  property  to  the  Laplace  Transform 
of  equation  2.3  yields 

•  *  \ 

Y (s)  =  G(s)  X (s)  2.6 


When  the  input  function  x(t)  is  a  unit  impulse  or  Dirac  delta  function 
6 ( t )  such  that 

/  6  (t)  dt  =  1 

—  00  V  ' 

2.7 


6(t) 


0  ,  t  i  0 


' 
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the  Laplace  Transform  of  x(t)  is  unity  and  so  equation  2.6  now  gives 


Y(s) 


G(s) 


Thus  the  Laplace  Transform  of  the  response  to  a  unit  impulse  is  the 
Laplace  Transform  of  the  weighting  function  g(t)  of  the  system  .This 
is  the  transfer  function  ,  G(s)  ,  of  the  system  which  is  only  defined 
for  a  linear  ,  time  invariant  and  noise-free  system  .  The  transfer 
function  completely  characterizes  the  input-output  relations  of  the 
system  . 

An  alternative  representation  of  the  transfer  function  is  the 
frequency  response  function  G(f)  given  by 


G  (f  ) 


2.3 


where  in  place  of  the  Laplace  Transform  we  are  using  the  one-sided 
Fourier  Transform  . 

The  frequency  response  function  is  complex  and  in  engineering 
contexts  its  magnitude  is  generally  plotted  as  the  dimensionless 
ratio  gain  in  dB.  .  However  ,  although  this  convention  is  often  used 


in  neurophysiology  (  e.g.  :  Poppele  and  Bowman,  1970  )  it  is 
inappropriate  when  the  input  and  output  signals  do  not  have  the  same 
dimensions  . 

When  the  input  x(t)  to  the  system  is  a  realization  of  a  stationary 


stochastic  process  one  can  define  the  autocovariance  function  p  (x) 

XX 

of  x(t)  as  (  Jenkins  and  Watts,  1963  ;  Lee,  1967  ) 


S 
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PXX(T>  =  ^  fiT  X(t)  x(  t  +T  )dt  2.9 

and  the  autocovariance  function  p  (x)  of  the  response  y(t)  as 

lira  1  .T  f  \  N  . 

pyy(T)  =  T^»  2T  /-T  y(t)  y(  C  +T  )dt 

and  the  forward  cross-covariance  p  (t)  from  the  input  to  the  output 

pxy(T)  =  T-”  It  /-T  X(t)  y(  C  +T  )dt  2,10 

These  covariance  functions  are  the  mathematical  expectation  of  the 
product  of  the  values  which  the  stationary  random  processes  assume 
at  instants  separated  by  x'.  i.e.  for  the  autocovariance  function  of 
x(t)  where  x(t)  is  a  realization  of  an  ergodic  process. 


PXX(T) 


E  {  x ( t )  x (  t  +x  )} 


They  are  related  to  the  corresponding  correlation  functions  y  (x)  , 

XX 

y  (x)  .  and  y  (x)  where  a  is  the  standard  deviation  by 
yy  xy 


y  (x) 

2  ,  x 

0  P 

XX 

X  XX 

y  (x) 

2  ,  N 

a  p  (x) 

yy 

y  yy 

y  (x) 

S3 

°  a  pv„( 

xy 

x  y  xy 

2.11 


The  autocovariances  and  autocorrelation  functions  are  even  functions 
of  t  ,  the  cross-covariance  and  cross-correlation  are  generally 
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odd  functions  of  t. 

From  the  Wiener-Khinchine  relation  (  Lee,  1967  ;  Thomas,  1969  ) 
one  can  obtain  the  power  spectra  of  the  input  and  the  output  and  the 
cross-power  spectrum  from  the  Fourier  Transforms  of  the  corresponding 
input  ,  output  and  cross-covariance  functions  . 

The  Fourier  Transform  F(f)  of  a  function  f(t)  is  defined  by  : 


F (f )  =  r  f(t)  e  j2lTft  dt  2.12 

—  00 

and  the  inverse  transform  by 

f(t)  =  \  T  F (f )  e^27Tft  dt  2.13 

277  -°° 


where  j  =  /-l. 

Thus  the  input  power  spectrum  S  (f)  ,  the  output  power  spectrum 


xx 


S  (f)  and  the  cross-power  spectrum  S  (f)  are  given  by: 

yy  xy 


S  (f) 
xx 


„oo  /  s  “1 2irf t  , 

/  P  (t)  e  dx 

-00  XX 


S  (f) 

yy 


00 

/  p 

yy 


— j  2irf  x 

(x)  e  j  dx 


2.14 


S  (f) 
xy 


00 

/  P 

-oo  'xy 


-j2TTfx  , 
(x)  e  dx 


Note  that  S  (f)  and  S  (f)  are  real  ,  even  functions  while  S  (f) 
xx  yy  y 

is  complex  and  odd  .  (f)  is  often  considered  in  terms  of  its 

real  part  ,  the  real-  valued  even  co-spectrum  C  (f)  and  its 
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imaginary  part  ,  the  real-valued  odd  quad-spectrum  (f)  where 


S  (f) 
xy 


C  (f)  -  jQ  (f) 

xy  J  xy 


2.15 


or  alternatively 

C  (f) 
xy 


=  S  (f)  +  S  (f) 

xy  yx 


Vf) 


=  j[  s  (f) 

xy 


S  (f)  ] 
yx 


where  S  (f)  is  the  backward  cross-spectrum  . 
yx  •  L 

These  spectra  are  two-sided  ,  that  is  they  are  defined  for 
negative  as  well  as  positive  frequencies  ,  Although  two-sided 
spectra  simplify  the  analysis  ,  in  most  applications  negative 
frequencies  do  not  have  a  physical  interpretation  ,  and  so  from  the 
two-sided  spectra  one  can  obtain  one-sided  spectra  defined  only 
for  positive  frequencies  by  folding  the  two-sided  spectrum  over 
the  zero  frequency  axis  so  that  the  one  sided  spectrum  is  given 


by 


R  (f) 

xx 


2S  (f) 
xx 


0 


0  <  f  < 


f  <  0 


2.16 


Since  the  cross-spectrum  is  odd  to  obtain  its  one-sided  spectrum 
the  two-sided  spectrum  has  to  be  folded  over  both  axes  . 

From  either  the  one  or  two-sided  spectra  one  can  obtain  the 
frequency  response  function  from  the  following  relations  : 


G(f) 


S  (f) 

_ 

S  (f) 
xx 


2.17 


s 


. 
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or 


G(f) 


2.18 


and 


I  G(f)  |2  Sxx(f) 


2.19 


Note  that  the  frequency  response  function  G(f)  can  be  interpreted  as 
a  complex  linear  regression  coefficient  for  obtaining  y(t)  from  x(t) 

(  Benignus,  1969a). 

2.3.  Application  to  linear  systems  with  intrinsic  noise  . 

If  the  response  y(t)  of  the  system  is  contaminated  by  additive 
intrinsic  noise  then  equation  2.3  does  not  hold  and  a  transfer  function 
may  not  be  defined  .  The  response  will  be  given  by 


/q  x(t)  h(  t  -t  )  dx  +  n(t) 


2.20 


y(t) 


where  n(t)  is  that  part  of  the  output  due  to  the  intrinsic  noise  and 
h(t)  is  the  weighting  function  of  the  linear  system  . 

In  this  case  a  describing  function  may  be  defined  whose  frequency 
response  function  H(f)  is  given  by 


H(f  ) 


/ 


oo 


0 


2.21 


Note  that  the  describing  function  does  not  completely  characterize 
the  input-output  relations  of  the  system  ,  it  provides  the  best  (in 


2  N  ‘ 


12 

the  sense  of  minimum  mean  squared  deviations  )  linear  model  for  the 
system  . 

If  the  input  signal  to  the  system  is  a  realization  of  a  broad¬ 
band  stochastic  process  one  can  define  the  coefficient  of  coherence 
(  Tick,  1963  ;  Enochson  and  Goodman,  1965  ;  Bendat  and  Pier sol,  1966  ; 
Benignus,  1969b  )  by 


y2(f) 


Alternative  definitions  of  coherence  have  been  proposed  .  Wiener 
(1930  )  defined  a  coefficient  of  coherency  as 


^  Sxy(f)  I  ,  S  (f).S  (f)  >  0 

- -  xx  yy 


S  (f).  S  (f) 
xx  yy 


o 


2.22 


S  (f).S  (f)  =  0 

xx  yy 


Y(f) 


S  (f) 
xy 


[  S  (f).S  (f)] 

xx  yy 


1/2 


and  this  is  a  complex  function  which  is  not  invariant  when  the  x(t) 
and  y(t)  are  subjected  to  linear  transformations  .  Foster  and 
Guinzy  (  1967  )  and  Hinich  and  Clay  (  1968  )  define  the  coefficient 
of  coherence  y(f)  as  the  positive  square  root  of  the  right  hand  side 
of  equation  2.22  ,  which  although  it  is  a  real  number  and  is  invariant 
when  the  x(t)  and  y(t)  are  subjected  to  linear  transformations  has 
been  described  (  Tukey,  1967  )  as ’an  unlikely  quantity'. 

The  definition  of  equation  2.22  has  been  used  because  of  its 
analogy  with  the  coefficient  of  determination  (  the  square  of  the 


correlation  coefficient  )  between  two  random  variables  .  The  coherence 


' 
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function.,  is  a  normalized  measure  of  the  linearity  of  the  relation 
between  x(t)  and  y(t)  and  since  the  spectral  matrix 


S(f) 


s 

xx 

S 

xy 


(f) 

(f) 


S 

yx 

S 

yy 


(f) 

(f) 


is  non-negative  definite  (  Koopmans,  1964  ) 


S 

xy 


(f) 


S 

xx 


(f).  s  (f) 


and  so 

0  <  y2 (f )  ^  1  2.23 

2 

y  (f)  represents  that  proportion  of  S  (f)  which  may  be  accounted  for 

2 

by  linear  regression  onto  S  (f)  ,  and  (  1  -  y  (f)  )  represents  that 

XX 

proportion  which  cannot  be  accounted  for  by  linear  regression  .  In 

this  case  this  is  the  proportion  of  S  (f)  due  to  S  (f)  ,  where 

yy  nn 

S  (f)  is  the  power  spectrum  of  the  noise  n(t)  . 
nn 

From  equation  2.20  it  can  be  shown  that  (  Amos  and  Koopmans,  1963  ) 


y2(f ) 


1 


1 


+ 


S 

nn 


(f) 


2.24 


S  (f)  | H(f ) | 2 

XX 

From  this  expression  Stein  and  French  (  1970  )  have  obtained  a  relation¬ 
ship  betx^een  the  coherence  function  and  the  information  transmission 
rate  of  the  system  .  Shannon  (  1948  )  derives  the  channel  capacity 


in  bits/second  of  a  noisy  channel  for  continuous  signals  as 


— 
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C  =  /  J?  log  (  1  +  S  (f)  )  df 

U  Z  ss 


2.25 


where  S  (f)  is  the  signal  spectrum  and  S  (f)  the  noise  spectrum  ,  both 
s  o  nn 

of  which  are  band-limited  to  W  .  From  equation  2.24  and  2.25  the 
rate  of  information  transmission  is  approximately  I  bits/second  where 


I  ^  -  /”  log2  (  1  -  y2(f )  )  df 


2.26 


An  important  property  of  the  coherence  function  is  its  invariance 
under  linear  operations  (  Koopmans,  1964  )  .  This  means  that  the 
physical  systems  used  to  measure  x(t)  and  y(t)  will  not  influence 
the  value  of  the  coherence  function  if  the  measuring  systems  are 
linear . 

2.4.  A.pplication  to  nonlinear  systems  with  intrinsic  noise  . 

For  a  zero  memory  nonlinearity,  that  is  ,  one  whose  output  at 
a  given  time  depends  only  on  the  input  at  the  same  time  ,  Booton 
(  1952  )  has  shown  that  for  a  given  input  function  x(t)  there  is 
an  equivalent  quasi-linear  element  K  (f)  such  that  the  mean  squared 
difference  between  the  response  of  the  nonlinear  element  and  the 
response  of  the  quasi-linear  element  is  minimized  .  Thus  analagous 
to  equation  2.20  the  response  of  the  nonlinear  system  is  related  to 
the  input  by 


y (t)  =*  /”  x(t)  k(t  -t  )  dx  +  tt(t) 


2.27 


15 

where  N(t)  is  due  to  the  nonlinear  terms  in  the  response  and  to  any 
intrinsic  noise  . 

In  this  case  the  describing  function  K(f)  may  be  defined  by 
equations  similar  to  those  of  section  2.3  ,  and  the  coherence 
function  provides  a  measure  of  that  part  of  S  (f)  which  is  linearly 

yy 

2 

dependent  on  S  (f)  and  (  1  -  y  (f)  )  provides  a  measure  of  that 

xx 

part  of  S  (f)  which  is  due  either  to  nonlinear  transformation  of 

yy 

S  (f)  or  to  intrinsic  noise  . 
xx 

Equation  2.27  is  not  a  general  representation  of  the  response 
of  a  nonlinearity..  A  more  general  representation  of  the  response 
of  a  nonlinear  system  is  : 

y (t)  =  $[  f  x(t)  k(t  -  t)  dr  +  n(t) ] 

where  n(t)  is  a  noise  term  and  $[•]  is  a  nonlinear  operator.  However, 
this  general  approach  is  not  very  tractable. 
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CHAPTER  3. 

ANALYTICAL  TECHNIQUES. 

3.1.  Hardware. 

All  experimental  data  (  analog  signals  and  spike  trains  )  was 
recorded  on  a  four  channel  FM  tape  recorder  (  Thermionic  T  3000  ) 
for  subsequent  analysis  on  a  Digital  Equipment  Corporation  LAB-8 
computer.  The  computer  configuration  used  consisted  of: 

1) .  an  8/1  central  processor  with  4  K  of  core  memory  and  an 
extended  arithmetic  element  KE  8/ I 

2) .  an  AX08  Laboratory  Peripheral  with  four  pulse  inputs,  a 
four  channel  multiplexed  analog-to-digital  convertor  (  ADC  ) ,  a 
crystal  and  an  RC  clock  and  a  storage  display  oscilloscope 

3) .  a  DF  32  magnetic  disc  unit  and  controller,  with  32  K  twelve 
bit  words  of  storage 

4) .  a  digital  input  buffer  capable  of  accepting  ten  pulse  train 
inputs  simultaneous  iy. 

The  principal  flow  paths  of  the  data  processing  system  are  given 
in  Figure  3.1.  Since  the  pulse  outputs  (  0  ->  +5V  )  from  the  pulse 
height  analyzer  units  (  Stein,  1968a)  are  incompatible  with  the  pulse 
input  sockets  of  the  AX08 ,  which  require  0  -*  -3V  ,  the  digital  input 
Buffer  was  used  for  all  spike  train  inputs  via  the  pulse  height  analyzer. 
A  further  advantage  of  using  the  digital  input  buffer  is  that  it  operates 
in  interupt  mode  ,  and  so  the  processor  can  perform  data  processing 
operations  while  accepting  data.  fhe  ADC  accepts  voltages  in  the  range 


s 
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Figure  3.1  .  Principal  information  flow-paths  of  data  processing 


system 
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+  1.024  V  full  scale.,  and  operates  by  successive  approximation  to 
nine  bits  .  Since  the  output  from  the  Thermionic  tape  recorder  is 
limited  to  +  1  V  the  analog  signals  could  be  fed  directly  to  the  ADC  . 

Data  output  from  the  system  was  either  via  the  storage  oscillo¬ 
scope  by  photography  or  via  the  teletype  . 

3.2.  Time  domain  processing  of  spike  trains. 

All  time  domain  processing  of  spike  trains  was  accomplished  using 
the  PULSE  and  FNEW-PULSE  programs  (  French,  1970  )  .  PULSE  is  a  general 
purpose  histogramming  program  which  accepts  two  pulse  trains  .  It  constructs 
a  histogram  with  a  variable  delay  (  0  to  0.4095  seconds  )  ,  a  variable 
number  N  of  bins  (  1  fo  128  )  of  bin  width  B.100  yseconds  ,  a  variable 
number  of  events/sweep  and  a  variable  number  of  sweeps.  These  terms  are 
defined  in  Figure  3.2  . 

FNEW-PULSE  provides  a  means  of  accessing  the  stored  histogram  by 
the  higher  level  interpretive  language  FOCAL  .  This  permits  fitting  the 
histogram  by  an  appropriate  curve. 

A  cycle  histogram  may  be  built  up  when  the  spike  train  is  the 
response  to  a  periodic  driving  function  .  Such  a  histogram  results  when 
the  trigger  pulse  is  synchronized  to  the  driving  function  and  the  product 
N.B  is  set  to  the  period  of  the  driving  function  ,  and  the  number  of 
events/sweep  is  large.  A  cycle  histogram  is  an  estimate  of  the  probability 
density  of  spikes  as  a  function  of  the  phase  of  the  periodic  driving 
function  . 

An  interspike  interval  histogram  may  be  built  up  when  the  number  of 
events/sweep  is  one  and  the  spike  train  is  serving  as  the  tiigger  pulse  . 

If  the  number  of  events/sweep  is  large  an  au to-correlogram  is  formed  . 


'  '  **  1 
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Figure  3.2  .  Diagram  of  example  histogram  defining  the  parameters 


of  section  3.2  . 
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3.3.  Methods  of  spectral  estimation. 

The  power  spectra  defined  by  equations  2.14..  are  obtained 
from  time  functions  f(t)  where  -°°  <t  <  +°°  .  In  practice  only  a  finite 
realization 


fT(t) 


f  (t)  0  <<  t  <  T 

0  T  <  0,  t  >  T 


3.1 


is  available  and  so  in  place  of  the  power  spectra  one  can  define  an 
estimate  of  the  power  spectrum  ,  the  sample  spectrum  (  Jenkins  and  Watts, 
1968  )  .  Thus  in  place  of  equation  2.14a  the  sample  spectrum  is  given  by 


S  (f) 

xx 


,T  ,  >.  -j2TTfx 

r-T  Pxx(t)  6  dT 


3.2 


where  p  (t)  is  a  sample  auto-covariance  function  given  by: 
xx 


Pxx(t)  =  yT  fT_T  x(t)  x(t  +t)  dt 


3.3 


An  alternative  definition  of  S  (f)  may  be  given  by: 

XX 


s  (f) 

XX 


1  fT  .  s  -i 2*7Tf t J  .T  ,  j2uft'  . 

=  -  /0  xT(t)  e  J  dt  /0  xT(t’)  e  dt' 


i  XT(f)  .  X*(f) 


Y  I  xT(f )  | 


3.4 


where  XT(f)  is  the  Fourier  Transform  of  xT(t)  and  X*(f)  is  its  complex 
conjugate. 

There  are  a  variety  of  methods  available  for  obtaining  the  sample 

spectrum  S  (f)  and  these  methods  may  be  considered  in  three  groups  : 
v  xx 


1)  Blacl<man-Tukey  methods 

2)  Direct  methods 

3)  Auto-regressive  methods. 
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alternative  classification  of  available  methods  is  given  in  Parzen  (  1968  ) 
The  Blackman-Tukey  ,  or  indirect  ,  methods  follow  from  equations  3.1 
and  3.2  (  Blackman  and  Tukey,  1958  ;  Jenkins  and  Watts,  1968  )  .  In  order 
to  obtain  an  estimate  of  the  spectrum  from  x^(t)  ,  the  auto— covariance 
function  p(x)  is  estimated  by  : 


P  (t) 
xx 


I  /Ht 

T  0 


0 


X  (t)  X  (t  +  | x | )dt ,0  <  I t I  <  T 


T  >  T 


3.5 


This  estimate  of  the  auto-covariance  function  is  only  considered  up  to  a 
a  lag  t  =  W  where  W  <<  T  .  Thus  the  auto-covariance  function  estimate 
has  been  multiplied  by  a  function  w  (t)  given  by  : 


wr(x) 


1 

0 


|t|  W  <<  T 
|  T  |  >  W 


3.6 


Alternatively  ,  the  auto-covariance  estimate  has  been  viewed  through  a 
rectangular  lag  window  of  width  2W.  Since  multiplication  in  the  time 
domain  is  equivalent  to  convolution  in  the  frequency  domain  this  means 

/v  * 

that  the  spectral  estimate  S  (f)  has  been  convolved  with  W  (f) 

xx.  r 

where  W^(f)  is  the  Fourier  Transform  of  w^(r)  given  by  : 


W  (f) 
r 


^  f 


3.7 
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and  so  the  spectral  estimate  is  given  from  equation  3.2  as  : 

Sxx(f)  =  /  1  Wr(g)  Sxx(f  -  g)  dg  3.8 

S  (f)  is  a  smoothed  spectral  estimate  since  the  effect  of  the  convolution 

XX 

A 

of  S  (f)  with  W  (f)  is  to  spread  the  component  at  a  frequency  f  over 

XX  IT 

(sin2irfW) /27TfW  for  -  <»  <f  <  «>  .This  spreading  of  the  component  at  f 
to  other  frequencies  is  called  leakage  from  the  frequency  f.  The  decay 

A 

of  the  leakage  away  from  the  frequency  f  is  shown  in  Figure  3.3  .  S  (f) 

XX 

is  an  unbiased  estimator  of  S  (f)  if  S  (f)  is  a  slowly  varying  function 

XX  XX 

of  f.  However,  if  the  spectrum  S  (f)  contains  large  fluctuations  there 

XX  ^ 

A 

will  be  considerable  leakage  and  hence  S  (f)  will  be  a  biased  estimator  . 

I  XX 

This  bias  may  be  reduced  by  quadratic  modification  ,  in  which  lag  windows 
other  than  w^(t)  are  used  .  A  desirable  property  of  these  lag  windows 
is  that  their  spectral  equivalents  should  have  lower  and  faster  decaying 
side  lobes  than  W^(f)  .  One  result  of  using  the  rectangular  lag  window 

A  A 

is  that  the  variance  of  S  (f)  is  2W/T  that  of  S  (f)  ;  since  W  <<  T 

xx  xx 

this  reduction  in  variance  is  considerable  ,  and  since  2W/T  tends  to 

A 

zero  as  T  tends  to  infinity  S  (f)  is  a  consistent  estimator  of  S  (f)  . 

xx  xx 

The  Blackman-Tukey  approach  has  been  extended  to  the  estimation  of 
cross-spectra  (  Jenkins  and  Watts,  1968  )  and  has  been  extensively 
modified  ,  notably  in  the  choice  of  the  optimal  lag  window  (  Akaike,  1962  ) 

•  x 

Discussion  of  the  statistical  properties  of  the  application  of  the 
Blackman— Tukey  approach  to  estimation  of  open— loop  frequency  response 
functions  have  been  given  (  Jenkins,  1963  ;  Akaike  and  Yamanouchi,  1962  , 

Akaike,  1964  )  . 

The  direct  method  of  spectral  estimation  follows  from  equation  3.4  , 
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Fourier  Transform  of  rectangular  lag  window  . 


t 


Figure  3.3  . 
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and  it  has  only  become  computationally  efficient  since  the  rediscovery 
of  the  Fast  Fourier  Transform  algorithm  by  Cooley  and  Tukey  (  1965  ), 
which  provides  an  economical  method  for  computing  the  Discrete  Fourier 
Transform  .  This  method  was  applied  for  the  estimation  of  all  spectra 

t 

and  so  will  be  described  in  greater  detail  in  subsequent  sections  of 
this  chapter  . 

Recently  Akaike  (  1969,  1970  )  and  Parzen  (  1968  )  have  argued  in 
favour  of  estimating  S  (f)  by  autoregressive  decomposition  .  Samples 

XX 

x (n)  of  the  time  series  are  assumed  to  have  the  form  : 
x  (n)  =  S  c  e(n-l) 

p=0  15 

where  c^  =1  and  the  e(n)fs  are  completely  uncorrelated  .  With  this 
assumption  x(n)  can  be  approximated  by 

x  *  P 

x(n)  =  e(n)  +  E  a  x(n-m) 

i  P»ra 

m=l  1 

and  if  x(n)  is  a  realization  of  a  weakly  stationary  process  this  approximation 

is  stable  and  can  be  made  as  accurate  as  desired  by  increasing  p  .  The 

magnitude  of  p  can  be  estimated  from  the  auto-covariance  function  estimate  . 

From  estimates  of  the  auto-regressive  coefficients  a  an  auto-regressive 

p,m 

estimator  of  the  normalized  sample  spectrum  may  be  obtained  .  It  has  been 
claimed  (  Gersch,  1970  )  that  this  method  gives  smoother  spectra  which 
are  more  readily  interpretable  than  spectral  estimates  obtained  by  the 
Blackman-  Tukey  method  or  the  direct  method  ,  and  that  this  method  is  more 
objective.  However,  this  method  has  as  yet  found  little  application  .  It 
appears  to  be  better  suited  for  the  estimation  of  the  spectra  of  noisy 
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signals  than  for  estimating  spectra  of  random  signals  and  spike  trains 
in  order  to  estimate  frequency  resonse  functions  or  coherence  functions  . 


3. A.  Sampling  a  continuous  signal  . 

The  computation  of  estimates  of  S  (f)  on  a  digital  computer  implies 

XX 

that  in  place  of  the  continuous  parameter  process  {  x(t)  ,  -°°<  t  <  00  } 


one  must  use  a  discrete  parameter  process  {  v(n)  :  n  =  0,  ±1,  ±2...  } 

which  may  be  formed  by  setting 


X(At)(n) 


x(nAt) 


choice  of  At  follows  from  the  Nyquist  sampling  theorem  (  Nyquist,  1928  ; 
Landau,  1967  )  .  If  x(t)  is  a  signal  of  finite  energy  bandlimited  to  B  hz. 

i.e. 


with 


x(t) 


I 


v  /  sin  27TB(t  -  k/2B  )  ■, 

-oo  \  1  2TTB(t  -  k/2B  )  ' 


3.9 


<  OO 


then  x ( t )  can  only  be  recovered  error-free  from  x^t)(n)  if 

1/At  >  2B 

Half  the  sampling  rate  ,  1/  2At  ,  is  called  the  Nyquist  or  folding 
frequency  f  .  If  x(t)  has  a  spectrum  S^Cf)  ,  then  it  can  be  shown 
(  Blackman  and  Tukey,  1958  )  that  SxxA(f)  ,  the  spectrum  of  x(At)(n)  Is 
given  by 


S  (f)  + 


xx 


[  S  (2kf  -f)  +  S  (2kf  +f)  ] 
L  xx  N-  xx  N 


3.10 
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Thus  if  there  is  any  power  in  S  (f)  for  f  >  f T  this  power  will 
alias  as  power  at  frequencies  f  <  f^  .  The  only  method  of  preventing 
such  aliasing  when  taking  regular  samples  is  to  choose  At  such  that 
>  all  where  S  ^(f^)  ^  0  .  This  may  be  simply  achieved  by 

filtering  x(t)  before  sampling  ,  and  choosing  an  appropriate  At.  When 
At  is  such  that  the  regular  samples  are  alias-free  ,  recovery  of  x(t) 
from  Is  stable  ,  that  is  ,  small  errors  such  as  quantizing 

errors  (  Gold  and  Rader,  1969  )  in  the  samples  give  small  errors 

in  the  recovered  x(t)  and  hence  in  estimates  of  S  (f)  . 


xx 

A  variety  of  methods  for  alias-free  sampling  of  a  stochastic  process 
have  been  proposed  in  which  the  sampling  times  t.  are  random  (  Beutler, 
1970  )  .  When  the  sampling  times  form  a  Poisson  process  ,  or  are 
obtained  by  random  deletions  from  a  regular  process  , alias-free  sampling 
can  be  achieved  at  mean  sampling  rates  less  than  the  Nyquist  rate  . 
However  such  techniques  are  not  readily  implementable  . 


3.5.  Sampling  a  spike  train. 

A  recorded  spike  train  is  a  ba  nd-limited  signal  of  finite  energy  as 
defined  by  equation  3.9  ,  but  since  the  duration  of  action  potentials 
is  of  the  order  of  milliseconds  regular  alias-free  samples  would  have  to 
have  an  excessively  high  sampling  rate  .  A  variety  of  methods  have  been 
proposed  for  sampling  a  spike  train  ,  and  these  have  been  criticized 
elsewhere  (  French  and  Holden,  1971a  )  together  with  a  full  description 
of  a  new  method  for  obtaining  alias— free  samples  of  a  spike  train  . 

In  the  method  used  to  obtain  alias-free  regular  samples  of  a  spike 
train  the  spike  train  is  treated  as  a  stochastic  point  process  •  Bartlett 
(  1963  )  has  defined  the  periodogram  (see  section  3.7  )  of  a  stochastic 
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point  process  as  the  Fourier  Transform  of  dn(t)  ,  where  dn(t)  is 

the  differential  of  the  cumulative  distribution  function  n(t)  .  n(t) 

of  any  point  process  is  a  staircase  function  ,  with  jump 

discontinuities  of  unity  height  at  the  occurance  times  of  the  events  . 

Thus  dn(t)  is  a  series  of  unit  impulses  or  Dirac  delta  functions 

at  the  occurance  times  t.  of  the  events  .  An  alternative  justification 

x 

for  treating  a  spike  train  as  a  series  of  unit  impulses  is  given 
by  Nelsen  (  1964  )  .  Thus  treatment  of  a  spike  train  as  a  stochastic 
point  process  is  identical  to  treatment  as  a  series  of  unit  impulses 
given  by 

i 

-  00 

y (t)  =  E  6(  t  -  t  )  3.11 

i=l 

in  order  to  estimate  the  periodogram  . 

The  spectrum  of  a  series  of  unit  impulses  extends  to  an 

i 

arbitrarily  high  frequency,  so  S  (f)  >  0  for  0  <  f  <  00  ,  and 

so  unless  the  impulse  train  is  filtered  aliasing  will  be  unavoidable  . 
Among  the  advantages  of  digital  filters  over  analog  filters  is  that  a 
digital  filter  can  have  an  impulse  response  x^hich  exists  for 
negative  times  .  Figure  3.4  shows  the  filter  characteristics  of  an 
ideal  ,  rectangular  low  pass  filter  ,  and  its  impulse  response  , 
which  is  given  by 

i 

g(t)  =  fN  [  sin(2rfNt)  ]/27TfNt 

3.12 

—  CO  <  t  <  00 

Such  a  filter  may  only  be  generated  digitally  . 


' 
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a 


\ 


Figure  3. A 


Ideal  ,  rectangular  low-pass  filter  and  its  impulse 


response  . 
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Thus  the  series  of  delta  functions  of  equation  3.11  may  be 
digitally  filtered  by  convolving  each  unit  impulse  with  g(t)  . 
Implementation  of  this  technique  may  be  considered  conceptually  as 
a  sequence  : 

i)*  f  y(t)  0  <  t  <  T 

yT(t)  =  |  . 

[  0  t  >  T 

ii) .  each  unit  impulse  at  t  is  convolved  with  g(t)  to  give  a 
function  y*  (t)  for  -<»  <:  t 


m)  . 


y’(t) 


iv)  . 


y’(t) 


0  <  t  <  T 


0 


t  >  T 


y^(nAt) 


Thus  the  result  is  equispaced  samples  of  the  digitally  filtered  train 
of  unit  impulses  ,  and  where  f  =  1/2 At  these  samples  are  alias- 

free  .  Details  of  practical  aspects  of  the  implementation  of  this 
technique  are  given  in  French  and  Holden  (1972  ). 


3.6.  The  Discrete  and  Fast  Fourier  Transforms  . 

The  sampling  methods  of  sections  3.4  and  3.5  generate  from 
xT(t)  and  y  (t)  a  series  of  N  numbers  ,  xT(At)(n>  and  yi(At)^  f°r 
n  =  0,  1,  ...N-l  .  The  Fourier  Transform  of  such  sets  of  numbers  may 
not  be  defined  .  However  ,  a  discrete  set  of  numbers  may  be 
associated  with  a  time  function 
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N-l 

M  =  E  x(nAt)  6(t  -  nAt) 
n-0 


whose  Fourier  Transform  is  given  by 


X’(f) 


=  r  x’(t) 
—00 


e 


“j2irf  t 


dt 


N-l 

Z 

n=0 


x (nAt) 


-j  2  Trf  At 
e  J 


3.13 


This  method  of  obtaining  the  Fourier  Transform  of  x(nAt)  provides 
a  means  of  defining  the  Discrete  Fourier  Transform  of  a  set  of  N 
numbers  x(n)  ,  n  =  0,  1,  _  N-l  as 


X(m) 


/  \  -i27rnm/N 
Z  x  (n)  e  J 

n=0 


3.14 


and  the  inverse  transform  by 


x  (n) 


1 

N 


N-l 

Z 

m=0 


X(m)  ej2™™/N 


3.15 


Thus  the  sequences  x(n)  and  X(m)  which  form  a  Discrete  Fourier  Transform 
pair  are  periodic  ,  with  a  period  N  .  This  may  be  visualized  by 
considering  x(n)  and  X(m)  to  be  defined  at  integer  points  on  a  circle 
of  circumference  N  (  Cooley, Lewis  and  Welch,  1969  ;  Gold  and  Rader, 

1969  )  .  Thus  the  Discrete  Fourier  Transform  is  not  simply  a  discrete 
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analog  of  the  Fourier  Transform  ,  but  has  some  properties  of  its 
own  . 

The  Fast  Fourier  Transform  ,  rediscovered  by  Cooley  and  Tukey 

(  1965  )  ,  is  an  algorithm  for  the  efficient  evaluation  of  the  Discrete 

Fourier  Transform  similar  to  that  described  by  Danielson  and  Lanczos 

2 

(  1942  )  .  The  evaluation  of  equation  3.14  requires  N  complex 

multiplications  and  divisions  ,  and  so  until  the  rediscovery  of  the 

Fast  Fourier  Transform  the  indirect  method  of  estimating  spectra  was 

computationally  more  efficient  than  the  direct  method  when  N  was 

large  .  There  are  two  kinds  of  Fast  Fourier  Transform  algorithm  , 

decimation  in  time  (  Cooley  and  Tukey  ,  1965  )  and  decimation  in 

frequency  (  Gentleman  and  Sande,  1966  )  .  Both  of  these  algorithms 

result  in  the  reduction  in  the  number  of  operations  needed  to  compute 

2 

an  N  point  transform  from  N  to  N  log£N  when  N  is  an  integral  power 
of  2  .  As  well  as  the  saving  in  computing  time  it  has  been  claimed 
that  the  Fast  Fourier  Transform  is  more  accurate  than  the  Discrete 
Fourier  Transform  since  the  smaller  number  of  operations  produce 
fewer  round-off  errors  (  Wienstein,  1969  )  . 

3.7.  Consistent  spectral  estimation  by  the  direct  method  . 

If  x (n)  :  n  =  0,  1,  ...,  N-l  are  N  alias-free  samples  of  x(t) 
taken  at  regular  intervals  At  from  a  record  of  length  T  =  NAt 
application  of  equation  3.14  gives  the  Fourier  coefficients  of 
x(nAt)  ,  X(mAf)  .  These  Fourier  coefficients  are  complex  and  so 


=  A  + 
m^ 


B 

m 


XCmAf) 
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where  the  index  Af  has  been  omitted  for  clarity  in  subsequent 
equations.  The  Fourier  coefficients  occur  at  equally  spaced  points 
with  a  frequency  resolution  of  1/NAt  . 

An  estimate  of  the  spectrum  ,  the  periodogram  I  (mAf)  is  given 

XX 


by 


I  (mAf) 
xx 


X (mAf)  X* (mAf) 

N 


X(mAf  ) 


N 


12  ? 

^  [  A  +  bJ  ] 

N  m  m 


3.16 


The  periodogram  is  purely  real  and  as  it  is  an  even  function  only 

[  N/2  ]  +1  values  are  distinct  .  It  is  not  a  good  estimate  of  the 

spectrum  S  (f)  since  it  is  biased  and  is  not  consistent  .  The  bias 
xx 

occurs  since  the  time  series  has  been  viewed  through  a  rectangular 
data  window  of  width  T  ,  and  so  the  coefficients  have  been  convolved 
with  T  [  sin  nfT  ]/  irfT  ,  and  so  the  periodogram  is  given  by 


I  (mAf) 
xx 


=  /  T 


sin  irf’T 


irfT  _ 


2 

S  (f-f')  df' 
xx 


3.17 


2 

where  T  [  sin  TTf’T/  tt f  tT  ]  is  the  Bartlett  spectral  window. 

In  the  case  where  x(t)  is  a  spike  train  which  does  not  have  a 
constant  interspike  interval  S  (0)  >>  S  (f  4-  0  )  ,  and  so  this 

XX  A.A 

convolution  will  produce  massive  leakage  from  the  BC  component  •  This 
DC  leakage  is  often  a  problem  and  so  is  always  eliminated  by  generating 
from  x (n)  :  n  -  0,  1,  ..  ,N-1  a  series  with  a  zero  mean  and  hence 


no  DC  component  . 


„J(  «  Ik  j.  »■  0  "  '  ’*  '• 

- 
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However  ,  there  will  still  be  bias  if  the  spectrum  S  (f)  is 

xx 

not  smooth  .  This  bias  due  to  leakage  may  be  reduced  by  linear 
modification  of  the  Fourier  coefficients  -  a  simple  and  widespread 
method  is  Hanning  ,  in  which  the  Fourier  coefficients  are  convolved 


-(1/ 4)A  3.18 

m+1 

and  similarly  for  the  B  .  Convolution  of  the  Fourier  coefficients 

m 

with  these  Hanning  weights  is  equivalent  to  viewing  the  time  samples 

/ 

through  a  cosine  bell  w^(n)  given  by 

w  (n)  =  1  [  1  -  cos  2irn  ] 

h  2 

N-l 


with  -  1/4  ,  +  1/2  and  -1/4  so  that 

Ao  "  1/2  1  A1  +  A0  1 

Am  =-(1/4)Vl  +(1/2)A„ 

an-i=  1/2  [  Vl  +  V2  1 


A  disadvantage  of  linear  modification  is  that  it  results  in  a 

decrease  in  frequency  resolution  with  only  a  small  decrease  in 

2  2 

variance  .  Quadratic  modification  ,  convolving  (A  +  B  )  with 

mm 

an  apparently  suitable  spectral  window  ,  has  little  effect  on 
leakage  (  Bingham,  Godfrey  and  Tukey  ,  1967  ;  Sloane,  1969  )  but 
does  permit  the  trading  of  a  decrease  in  frequency  resolution  for 
a  decrease  in  variance  . 

Obtaining  a  zero-mean  series  and  Hanning  can  reduce  bias  due  to 
leakage  ,  but  the  linearly  modified  periodogram  is  not  a  consistent 
estimator  of  the  spectrum  .  If  x(n)  are  samples  from  a  Normal  process  , 


N 
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since  the  Fourier  Transform  is  a  linear  transformation  A  and  B 

m  m 

will  be  Normally  distributed  and  so  I  (mAf)  will  have  a  chi- 

xx 

squared  distribution  with  two  degrees  of  freedom  .  If  x(n)  are 


samples  from  a  process  which  does  not  have  a  Normal  distribution  but 


N  is  large  then  from  the  Central  Limit  Theorem  A  and  B  will  have 

m  m 

approximately  Normal  distributions  and  so  I  (mAf)  will  have 

xx 

approximately  a  chi-squared  distribution  with  two  degrees  of  freedom  . 


The  expected  value  of  I  (mAf)  is 

xx 


E  {  I  (mAf)  } 
xx 


=  /  T r sinnf 'T,  S  (f-f')  df ' 

“°°  [ - FwH  xx 

i  -n-f'T 


3.19 


and  the  variance  is  approximately  S  (f)  . 

XX 


If 


W(m) 


N-l 

E  w,  (n)  e 
n=0 


-j  2irnm/N 


3.20 


where  w,  (n)  is  the  Hanning  window  then  the  expected  value  of  the 
h 

linearly  modified  periodogram  I  (mAf)  is 

XX 


E  {  I  (mAf)} 
xx 


2  E{  X(m)  }  E  W (ra) 2 


N 


3.21 


and  the  variance 


Var  {  I  (mAf)} 
xx 


4  E2{  X(m)2}  [  W (m) 2  ]2 

N2 


3.22 


Var  (I  (mAf ) }  =  E2  {  Ixx(mAf)} 


xx 


hence 


" '  .  u 
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as  in  the  case  of  the  unmodified  periodogram  (  Jenkins  and  Watts,  1968  ) 
The  variance  of  the  modified  periodogram  I  (mAf)  may  be  reduced 

XX 

by  averaging  M  periodogram  estimates  from  non-overlapping  records 

A 

of  length  T  .  The  variance  of  the  averaged  modified  periodogram  I  (mAf) 

XX 

is  now  approximately 


s  4(f) 


and  the  averaged  ,  modified  periodogram  has  approximately  a  chi- 
squared  distribution  with  2M  degrees  of  freedom  .  A  further  decrease 
in  variance  to  55%  can  be  achieved  by  having  the  record  segments 

A 

overlap  by  T/2  (  Welch,  1967  )  .  If  the  resolution  of  I  (mAf)  is 
unnecessarily  fine  a  further  reduction  in  variance  can  be  achieved 
by  averaging  adjacent  frequency  components  (  Jones,  1965  )  . 

When  k  is  the  number  of  effective  degrees  of  freedom  of  the 
averaged  estimate  of  the  modified  periodogram  Jenkins  and  Watts  have 
shown  that  100 (  1  -  a)%  confidence  limits  for  the  estimate  are  given 
by 


k 


I  (mAf) 


zfe(  1  -  a/2) 


S 

xx 


(mAf) 


k  I  (mAf) 
xx 

*k(a/2) 


where  z  (  1  -  a/2)  and  z  (a/2)  are  the  a/2  and  1  -  a/2  points  on 
k  k 

the  cumulative  chi— squared  distribution  with  k  degrees  of  freedom  as 

i 

given  in  Bendat  and  Piersol  (  1966  )  . 


s 

1 
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3.8.  Consistent:  cross-spectral  estimation  by  the  direct  method  . 
The  sample  cross-spectrum  is  defined  by 


S  (f) 
xy 


-  X*(f)  YT(f) 


3.24 


”  [Ax(f)  -  jBx(f)  ]  [  Ay(f)  +  jBy(f)  ] 


where  [  A  (f)  +  jB  (f)  ]  and  [  A  (f)  +  jB  (f)  ]  are  the  Fourier 

x  x  y  y 

Transforms  of  x^(t)  and  y  (t)  .  Estimation  of  the  cross-spectrum  by 

A 

the  cross-periodogram  I  (mAf)  is  analogous  to  the  estimation  of 

A 

S  (f)  by  I  (mAf)  :  the  record  is  broken  into  subrecords  of  length 

XX  XX 

T  ,  and  for  each  subrecord  ~  X*(mAf)  Y^(mAf)  is  evaluated  from 
linearly  modified  Fourier  coefficients  ,  and  periodograms  from  different 

A 

subrecords  are  averaged  to  form  I  (mAf)  . 

xy 

A 

I  (mAf)  can  be  considered  either  in  terms  of  its  magnitude  and 
xy 


phase 


where 


and 


I  (mAf) 
xy 


M  (mAf) 
xy 


P  (mAf) 
xy 


*  ,  A,s  jP  (mAf) 

M  (mAf)  eJ  xy 
xy 


I  I 


3.25 


tan 


-1 


IM  I  (mAf) 

xy  


RE  I  (mAf) 
xy 


or  in  terms  of  its  real  and  imaginary  parts  ,  which  are  estimates  of 
the  co-  and  quad-spectrum  . 

If  I  (mAf)  is  the  average  of  M  cross-periodogram  estimates  which 
xy 

have  not  been  subject  to  linear  modification  then  the  variance 


' 


37 

A 

of  (mAf)  is  given  by  Hinich  and  Clay  (  1968  )  as 


A 

Var{  M  (mAf)} 
xy 

’  IM  Sxx(f)  Syy(f)  [  1  ^2(f)  1 

and 

A 

Var{  P  (mAf)} 
xy 

1  [  1  -  1  ] 

2M  1  2  1  J 

Y  (f) 

Thus  the  variance  of 

the  cross-spectrum  estimate  is  dependent  on  the 

"\ 

coherence  function. 

3.9.  Consistent  estimation  of  the  coherence  function  . 

If  the  coherence  function  defined  by  equation  2.22  is  estimated 
from  the  linearly  modified  but  unaveraged  periodograms  by 


y2(f) 

I  Vf)  I2 

xy  3.26 

i  (f)  i  (f) 

xx  yy 

A  A 

it  will  be  identical  with  unity  everywhere  except  when  I  (f).I  (f) 

xx  yy 

is  zero  since  if  (  a  +  jb  )  ,  (  a’+  jb')  are  the  Fourier  coefficients 


of  xT(t)  and  yT(t) 

1  i  (£)  I  2 

1  xy  ' 

|  (  a  -  jb  )  (  a'  +  jb'  )  |  2 

T 

i  (f).  i  (f) 

xx  yy 

(  a  -  ib  )(  a  +  jb  )  (  a’-jb')(  a'+jb'  ) 

T  T 

=  1  . 

2 

Thus  equation  3.26  is  a  useless  estimator  of  y  (f)  ,  and  it  is 


•  t  .  V 


■ 

• 

' 

* 
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analogous  to  estimating  the  correlation  coefficient  from  a  single 
observation  pair  . 

A  widely  used  estimator  of  the  coherence  function  is 


Y  2  (f  ) 


(f) 


I 

xx 


(f) 


I  (f) 

yy 


3.27 


but  this  is  a  biased  estimator  .  Enochson  and  Goodman  (  1965  )  have 
shown  that  estimates  of  the  coherence  function  are  positively  biased 
and  that  this  bias  is  a  function  of  the  effective  degrees  of  freedom 
of  the  estimate  and  the  true  coherence  ,  and  have  given  curves  of  a 
bias  correction  factor  when  the  true  coherence  is  greater  than  0.4  . 
Benignus  (  1969b)  used  Monte  Carlo  methods  to  investigate  this  bias 
and  produced  an  empirical  family  of  curves  relating  the  bias  to  the 
effective  degrees  of  freedom  of  the  estimate  and  the  true  coherence 
for  all  values  of  coherence  between  zero  and  one  .  These  curves 
were  fitted  by  a  least  squares  procedure  and  the  estimate  of  bias 
B  [y2(f)  ]  is  given  by 

B  [  y 2 ( f )  1  =  l  (  1  -  Y2(£)  )  3.28 

and  so  that  the  corrected  estimate  of  coherence  is  given  by 

~9  ~2  ~  ~2 

y2(f)  =  y  (f)  -  B  [  y  (f)  ] 

A 

~2  2 

The  residual  error  between  y  (f)  and  y  (f)  was  less  than  0.01  for 
all  estimates  with  k  >  8  . 


■ 


. 
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Tick  (  1967  )  has  proposed  an  alternative  family  of  estimators 
2 

for  y  (f)  .  In  place  of  just  averaging  the  linearly  modified 

per iodograms  he  proposes  that  coherence  function  estimates  obtained 

from  averaged  modified  periodograms  should  themselves  be  averaged  . 

This  estimate  is  more  biased  when  the  true  phase  between  the  two 

signals  is  changing  slowly  ,  but  is  less  biased  when  the  true  phase 

is  changing  rapidly  .  In  Monte  Carlo  studies  ,  where  the  true 

coherence  is  known  ,  the  number  of  periodograms  averaged  to  be  used 

to  form  a  coherence  function  estimate  ,  and  the  number  of  coherence 

function  estimates  to  be  averaged  can  be  optimized  to  form  an  accurate 

estimate  .  However  ,  when  the  true  coherence  is  unknown  this  method 

is  not  very  useful  but  does  suggest  that  thepe  is  no  optimal  method 

for  estimating  an  unknown  coherence  function  . 

Benignus  (  1969a)  discussed  the  estimation  of  the  coherence 

function  by  use  of  the  Fast  Fourier  Transform  ,  and  favours  an 

estimator  produced  by  both  averaging  a  set  of  periodograms  and 

averaging  over  frequency  segments  within  the  periodograms  .  When 

averaging  over  frequency  segments  the  adjacent  components  are  not 

independent  if  the  coefficients  have  been  linearly  modified  ,  and  so 

care  must  be  taken  in  evaluating  the  effective  degrees  of  freedom  k  . 

Goodman  (  1957  )  has  derived  a  distribution  function  for  the 

positive  square  root  of  the  coherence  function  between  two  Gaussian 

processes  ,  and  this  has  been  tabulated  by  Amos  and  Koopmans  (  1963  )  . 

Using  Monte  Carlo  methods  Foster  and  Guinzy  (  1967  )  have  shown  that 

this  distribution  is  applicable  to  non-Gaussian  processes  if  the 

spectra  S  (f)  ,  S  (f)  and  S  (f)  are  reasonably  smooth  ,  and  by 

xx  yy 

using  a  Bayesian  approach  (  Good,  1965  )  have  obtained  a  maximum 
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likelihood  estimator  for  the  coherence  function  :  however  ,  when 

/V 

~2 

k  is  large  this  estimator  does  not  differ  significantly  from  y  (f)  . 

Enochson  and  Goodman  (  1965  )  have  shown  that  the  sampling 
distribution  of  the  Fisher  Z-transform  (  Kaplan  ,  1962  )  of  the 
positive  square  root  of  the  coherence  function  estimate  is  approximately 
Gaussian  with  a  variance  given  by 


A 

Var  {  Z  /y2 (f )  } 


1 

2k 


/V 

when  0.4  <  y2(f)  <  0.95  and  k  >20. 

Hinich  and  Clay  (  1968  )  have  given  the  variance  of 


/y2(f)  as 


Var  {  /y2(f)  }  a,  fk  t  1  -  A2(f)  1  2 

and  so  the  variance  of  the  coherence  function  estimate  depends  on  the 
value  of  the  true  coherence  . 

A 

~ 2 

Approximate  confidence  limits  for  y  (f)  are  given  by 


U(f) 


m{  1  -  a/2  } 

/  k 


3.29 


wher  m  {1  -  a/2)  is  the  1  -  a/2  point  on  the  cumulative  unit  Normal 
distribution  function  and 

A 

a 

U(f)  =  tanh  y  (f ) 


which  is  assumed  (  by  appealing 


to  the  Central  Limit  theorem  )  to  have 
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a  Normal  distribution  (  Bendat  and  Piersol,  1966  )  . 


3.10.  Consistent  estimation  of  the  frequency  response  function  . 

The  frequency  response  function  G(f)  defined  by  equations  2.16 
2.18  may  be  estimated  from 


Gx(f)  = 


or 


C2(f)  = 


I 

xy 


(f) 


I 

XX 


(f) 


I 

yy 


i  (f) 

yx 


Note  that 


GL(f) 

V) 


~2 

Y  (f) 


3.30 


3.31 


and  so  these  two  estimates  will  only  be  the  same  if  the  coherence 

A 

function  is  one  .  Bendat  (  1962  )  has  argued  that  G^f)  should  be 
used  to  estimate  G(f)  when  the  output  is  contaminated  with  noise  , 


and  G  (f)  should  be  used  when  the  input  measurement  system  introduces 
noise  .  G^(f)  was  used  to  estimate  G(f)  since  the  coherence  function 

was  never  one. 

A 

G(f)  is  complex  ,  and  confidence  limits  for  its  magnitude!  G(f)  | 
and  phase  /G(f)  are  given  by  Jenkins  and  Watts  (  1968  )  as 
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where  f 2 ,1^.-2  ^  a  point  on  the  cumulative  Fisher’s  F 

distribution  with  2  degrees  of  freedom  and  k-2  degrees  of  freedom 

A  A 

The  variance  of  |  G(f)  |  and  /G(f)  are  given  by 
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and  so  for  low  values  of  coherence  the  variance  is  high 


3.11.  Software  for  spectral  estimation  . 

The  input  signals  were  sampled  using  SAP-1  (II)  (  French  and 
Holden*  1972  )  and  the  Fourier  coefficients  obtained  by  SAP-3  which 
utilizes  a  time  decimation  Fast  Fourier  Transform  algorithm  (  French 
and  Holden,  1971b  )  .  Hanned  spectral  estimates  were  obtained  from 
the  Fourier  coefficients  and  averaged  by  SAP-9  (  French,  1971  )  and 
estimates  of  the  spectra  ,  frequency  response  functions  ,  coherence 
function  and  information  transmission  rate  and  their  confidence 


« 


43 


limits  were  generated  by  FOCAL  programs  accessing  the  stored  spectral 
estimates  by  SAP-4  (  French  and  Holden,  1971b  )  . 

A  flow  chart  showing  the  principal  computational  routes  is 

i*  . 

given  in  Figure  3.5  .  A  full  discussion  of  computational  aspects  of 
spectral  estimation  by  the  Fast  Fourier  Transform  are  given  in 
Cooley,  Lewis  and  Welch  (  1967  )  and  Brumbach  (1963  )  . 
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Figure  3.5  .  Flow  chart  of  principal  computation  routes  of  spectral 


analysis  package 
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CHAPTER  4  . 

THE  COCKROACH  TACTILE  SPINE  . 

4.1.  Introduction  . 

Pringle  and  Wilson  (  1952  )  utilized  linear  systems  theory  to 
obtain  a  quantitative  model  for  the  input-output  relations  of  the 
large  tactile  spine  on  the  dorsal  surface  of  the  femur  of  the  meta- 
thoracic  leg  of  the  cockroach  ,  Periplaneta  americana  .  Sinusoidal 
and  rectangular  wave  driving  functions  were  applied  to  the  receptor  , 
and  a  ’transfer  function’  relating  the  input  tension  to  the  rate 
of  sensory  spike  discharge  was  obtained  .  This  receptor  has  been 
studied  using  a  wide  variety  of  techniques  (Chapman  and  Smith  ,  1963  ; 
Crowe,  1967  ;  Pearson  and  Holden,  1970  ;  Holden  and  French,  1971a  , 
1971b  ;  French,  Holden  and  Stein,  1972  )  and  thus  it  provides  a  very 
convenient  preparation  for  the  development  and  comparison  of  different 
techniques  of  analysis  . 

Pringle  and  Wilson  took  the  rate  of  sensory  spikes  as  the  output 
from  the  receptor  ,  and  fitted  the  response  to  a  square  wave  tension 
driving  function  by  the  linear  sum  of  three  decaying  exponential  . 
functions  .  This  would  give  a  frequency  response  function  which  would 
have  in  phase  modulation  at  low  frequencies  ,  and  an  increasing 
sensitivity  and  phase  lead  passing  through  local  maxima  with  increasing 
frequency  .  However  ,  a  range  of  sinusoidal  driving  functions  adequate 
to  test  this  model  was  not  available  . 
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Chapman  and  Smith  (  1963  )  re-examined  the  preparation  using 
similar  driving  functions  ,  and  fitted  the  rectangular  wave  response 
by  a  power  function  decay  .  The  frequency  response  function  predicted 
from  this  response  has  a  constant  phase  lead  and  the  sensitivity 
increases  as  a  power  of  the  forcing  frequency  .  This  frequency 
response  function  was  supported  by  the  use  of  sinusoidal  driving 
functions  with  frequencies  from  0.1  to  lOhz.  . 

Crowe  (1967  )  compared  the  response  predicted  by  both  models  in 
response  to  ramp  driving  functions  with  the  experimentally  observed 
response  ,  and  his  results  tended  to  support  Pringle  and  Wilson’s 
model  . 

This  lack  of  agreement  suggests  that  the  frequency  response 
function  of  this  receptor  should  be  estimated  using  a  number  of 
techniques  and  the  results  compared  .  If  the  system  is  non-linear 
an  estimate  of  the  frequency  response  function  based  on  one  kind  of 
driving  function  may  not  necessarily  predict  the  response  to  other 
driving  functions  . 

4.2.  Methods  . 

All  experiments  were  performed  on  the  single  large  tactile  spine 
located  on  the  dorsal  surface  of  the  end  of  the  femur  of  either 
metathoracic  leg  of  adult  cockroaches  ,  Periplaneta  americana  .  The  leg 
was  cut  through  at  the  coxa  ,  and  the  cut  surface  covered  with  Vaseline 
£o  prevent  drying  .  The  leg  was  mounted  on  the  armature  of  a  Pye  Ling 
vibration  generator  (  three  ohm  input  impedance  ,  model  V47  )  by  two 
insect  pins  on  a  piece  of  cork  >  and  the  pins  served  as  recording 
electrodes  .  This  simple  recording  method  is  effective  because  of 
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the  high  resistance  of  the  narrow  femur  and  since  each  spine  is 
innervated  by  a  single  axon  (  Chapman  and  Nichols,  1969  )  . 

The  vibration  generator  was  driven  by  the  output  from  a  force 
servo  control  unit  via  a  power  amplifier  (  for  impedance  matching  )  . 

The  feedback  transducer  to  the  servo  was  the  differential  output 
from  two  Endevco  Laboratories  Pixie  8101  force  transducers  .  The  diff¬ 
erential  output  was  taken  in  order  to  compensate  for  temperature 
dependent  DC  shifts  .  The  transducers  respond  linearly  up  to  40  gms.-wt.  , 
however  ,  the  addition  of  a  lever  extension  limited  this  range  to  5 
gms.-wt.  .  The  extension  of  the  force  transducer  was  held  against  the 
spine  .  The  command  signal  to  the  servo  could  be  : 

i) .  the  output  from  a  Hewlett-Packard  3300A  function  generator 
which  had  passed  through  a  calibrated  amplitude  and  offset  control  unit 

ii) .  the  output  from  a  band-limited  white  noise  generator 

iii) .  any  linear  sum  of  (i)  and  (ii)  above  , 

The  amplified  action  potentials  ,  the  output  from  the  force 
transducer  and  reference  markers  were  recorded  on  a  Thermionic  Products 
T3000  FM  tape  recorder  for  subsequent  processing  on  the  LAB-8  computer  . 

Estimates  of  the  frequency  response  function  were  obtained  by 
subjecting  the  tension  servo  which  stimulated  the  receptor  to 

a)  deterministic  driving  functions  (sinusoids  and  rectangular  waves  ) 

i 

b)  sinusoidal  driving  functions  added  to  a  band-limited  white 
noise  auxiliary  signal 

c)  band-limited  white  noise  . 

Cycle  histograms  of  the  response  to  deterministic  sinusoidal 


48 


driving  functions  were  fitted  using  a  least  mean  squares  procedure 
(  Sokoln  ikof f  and  Redlieffer,  1958  )  .  When  band-limited  white  noise 
was  the  driving  function  the  input  ,  output  and  cross-spectra  were 
estimated  .  The  algorithms  and  programming  methods  used  have  been 
described  elsewhere  (  French,  1970  ;  French  and  Holden,  1971a  ;  1971b  ; 
1972  )  . 

i  ' 

4.3.  Results  . 

The  cycle  histograms  in  response  to  a  0.05  hz .  rectangular  wave 
driving  function  of  two  different  amplitudes  are  shown  in  Figure  4.1. 
These  histograms  could  not  be  fitted  by  a  single  decaying  exponential 
function  .  The  best-fitting  (  in  the  sense  of  least  mean  squared 
deviation  )  exponential  decayed  slower  than  the  histogram  in  the 
first  second  and  faster  than  the  histogram  in  the  last  few  seconds  . 

The  histogram  could  be  fitted  by  the  sum  of  several  exponentials 
using  a  method  of  successive  approximations  ,  however  this  was  not 
reproducible  .  A  log-log  plot  of  the  histogram  (  insert  on  Figure  4.1  ) 
could  be  fitted  by  a  straight  line  ,  with  a  slope  of  0.82  +  0.03 
and  a- linear  correlation  coefficient  of  0.80  +  0.01  for  fourteen 
different  spines  .  This  implies  that  the  histogram  may  be  fitted 
by  a  power  function  decay  of  the  form 

y ( t )  =  b  t  4.1 

where  y(t)  is  the  number  of  spikes/second  and  k  is  a  constant  given 
by  the  slope  of  the  log-log  plot  .  Although  minimizing  mean  squared 
error  on  a  log-log  plot  does  not  necessarily  minimize  mean  squared 
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Figure  4.1 
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Cycle  histograms  of  the  response  of  a  tactile  spine 

i 

to  a  0.05  hz.  rectangular  driving  function  of  two 
different  amplitudes  .  Inserts  show  log-log  plots  . 
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error  on  the  linear  plot  of  the  histogram  the  computed  curve  fits  the 
histogram  well  . 

This  step  response  is  similar  to  that  given  by  Chapman  and  Smith 
(  1963  )  for  the  decay  of  instantaneous  spike  rate  .  The  describing 
function  predicted  from  such  a  step  response  is  given  by 

H(f)  =  br(l-k)  fk  e  ^k7r/2  4 


where  T  is  the  gamma  function  defined  by 


r  (x) 


/ 


o 


u 


x  -  1 


-u 


du 


x  >  0 


4.3 


This  describing  function  has  a  constant  phase  advance  of  k-rr / 2  radians 

til 

and  the  sensitivity  of  the  receptor  increasing  as  the  k  power  of 
the  driving  frequency  . 

The  cycle  histogram  of  the  response  to  a  one  hz.  sinusoidal 
driving  function  is  not  a  sinusoid  -  Figure  4.2  shows  it  to  be  a 
clipped  ,  distorted  sinusoid  .  Hence  the  relation  between  the  tension 
driving  function  and  the  cycle  histogram  is  nonlinear  .  This  may  be 
further  illustrated  by  plotting  the  Fourier  coefficients  of  the  cycle 
histogram  .  A  linear  system  would  only  have  a  significant  component 
at  the  driving  frequency  (  assuming  that  there  was  no  output  in 
the  absence  of  an  input  )  .  In  fact  there  is  considerable  harmonic 
distortion  ,  for  where  harmonic  distortion  is  defined  by 
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Figure  4.2 
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Cycle  histogram  of  the  response  of  a  tactile  spine 
to  a  1  hz.  sinusoidal  driving  function  and  a  plot 
of  the  harmonic  Fourier  coefficients  of  the  histogram  . 
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where  a  and  b  are  the  real  and  imaginary  Fourier  coefficients  at 

n  n 


a  frequency  n  times  the  driving  frequency  ,  the  harmonic  distortion 
is  always  greater  than  one  .  The  presence  of  appreciable  first 


(  n  -  2  )  and  third  (  n  =  4  )  harmonics  in  Figure  4.2  is  compatible 


with  a  major  nonlinearity  being  that  of  rectification  . 

However  ,  the  amplitude  of  the  component  at  the  fundemental  ,  and 

the  amplitude  of  the  harmonics  ,  increase  linearly  with  the  amplitude 

of  the  sinusoidal  driving  function  at  this  frequency  .  This  linear 

relationship  ,  shown  in  Figure  4.3  ,  is  compatible  with  Crowe’s  (  1967  ) 

observation  that  the  maximum  instantaneous  rate  of  action  potentials 

was  linearly  related  to  the  peak  to  peak  tension  of  the  driving 

function  at  several  low  frequencies  of  sinusoidal  driving  functions  . 

Figure  4. ^  is  similar  to  that  observed  with  an  electronic  neuronal 

analog  (  French  and  Stein,  1970b  )  in  the  range  of  input  amplitudes 

which  produce  rectification  (  Stein,  French  and  Holden,  1972  )  .  This 

linear  relationship  provides  a  justification  for  estimating  the 

frequency  response  function  by  fitting  the  cycle  histograms  in  response 

to  sinusoidal  driving  functions  by  the  best-fitting  sine  wave  with  a 

frequency  equal  to  the  frequency  of  the  driving  function  .  The 

2  2 

amplitude  of  the  best— fitting  sine  is  J  [  a^  +  b^  ]  ,  a  measure  of 
the  power  at  the  driving  frequency  ,  not  the  peak  to  peak  amplitude 


which  has  been  widely  used  (  Borsellino,  Poppele  and  Terzuolo  ,  1965  ; 
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Figure  A. 3 
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Amplitude  of  the  fundemental  (  filled  circles  )  and 
the  first  harmonic  (  open  circles  )  of  cycle  histograms 
of  the  response  of  the  tactile  spine  to  sinusoidal 
driving  functions  of  increasing  amplitude  . 


54 


Melvill  Jones  and  Milsum,  1970  )  .  A  plot  of  the  logarithm  of  the 
amplitude  ,  and  the  phase  of  the  fitted  sine  curve  both  as  a  function 
of  the  logarithm  of  the  driving  frequency  form  a  Bode  plot  Xizhich 
characterizes  the  frequency  response  function  of  the  receptor  (  D’Azzo 
and  Houpis  ,  1966  )  . 

The  cycle  histograms  in  response  to  sinusoidal  driving  functions 
at  higher  frequencies  exhibit  phase  locking  -  the  action  potentials 
are  entrained  to  specific  phases  of  the  driving  function  .  Figure  4.4 
shows  the  cycle  histograms  in  response  to  sinusoidal  driving  functions 
at  5hz .  and  20  hz.  .  At  5hz.  there  is  a  strong  tendency  for  there  to 
be  five  spikes/cycle  occuring  at  phases  fixed  relative  to  the  driving 
function  ,  and  at  20  hz.  there  is  one-to-one  phase  locking  .  Phase 
locking  has  been  observed  in  a  variety  of  sensory  neurones  (  Kiang  et 
al.  ,  1965  ;  Matthews  and  Stein,  1969a)  and  in  simple  neuronal  models 
(  Rescigno  et  al.  ,  1969  ;  Stein,  French  and  Holden,  1972  )  .  For 
frequencies  where  one-to-one  entrainment  occurs  the  amplitude  of  the 
fitted  sine  curve  will  be  twice  the  driving  frequency  . 

A  Bode  plot  of  the  frequency  response  function  obtained  by 
using  sinusoidal  driving  functions  is  given  in  Figure  4.5  ,  and  is 
compared  with  the  frequency  response  function  predicted  from  the  step 
response  of  Figure  4.1  by  equation  4.2  .  The  agreement  between  the 
predicted  curve  (  solid  line  )  and  the  data  points  is  good  up  to 
lOhz .  ,  above  which  the  slope  of  the  sensitivity  curve  obtained  by 
the  sine-fit  method  is  less  than  that  predicted  from  the  step 
response  .  This  is  to  be  expected  as  phase  locking  occurs  at  frequencies 
above  10  hz.  . 

Thus  the  use  of  deterministic  input  signals  shows  two  kinds  of 


v 


55 


. 


.1X1 


lift! 


rdul 


Figure  4.4  .  Cycle  histograms  of  the  response  of  the  tactile  spine 

to  5  hz .  and  20  hz.  sinusoidal  driving  functions  showing 


phase  locking  . 
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Figure  4.5 


Bode  plot  of  the  frequency  response  function  of  the 
tactile  spine  obtained  by  using  sinusoidal  driving 


functions  . 
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nonlinearity  when  the  cycle  histogram  is  taken  as  the  output  - 
rectification  at  all  frequencies  and  phase  locking  at  high  frequencies  . 
From  the  cycle  histogram  in  response  to  sinusoidal  driving  functions 
the  harmonic  distortion  can  be  measured  ,  however  this  is  not  a  good 
estimate  of  the  linearity  of  the  system  as  no  consideration  is  taken 
of  nonharmonic  distortion  .  Bayly  (  1968  )  has  shown  that  even  in 
the  case  of  a  perfect  or  linear  voltage  to  rate  convertor  there  is 
appreciable  nonharmonic  distortion  even  with  small  input  signals  . 

Such  nonharmonic  terms  may  not  be  recovered  from  the  cycle  histogram  , 
as  for  a  driving  function  with  a  period  T  the  cycle  histogram  is 
viewed  through  a  spectral  window  of  shape  sin  (TrfT)/irfT  as  illustrated 
in  Figure  3.3  when  T  =  2W  .  This  leakage  from  the  coefficients  will 
bias  estimates  of  coefficients  at  all  other  frequencies  except  the 
harmonics  .  Estimates  of  the  coefficients  at  the  fundemental  and  the 
harmonics  will  not  be  biased  by  leakage  from  nonharmonic  coefficients 
as  building  up  a  cycle  histogram  has  the  effect  of  averaging  out  any 
coefficients  which  do  not  have  a  constant  phase  relation  with  the 
driving  function  . 

Rectification  occurs  in  the  cycle  histogram  since  there  is  no 
carrier  rate  of  spikes  which  can  be  modulated  .  Rectification  can 
also  occur  when  there  is  a  carrier  rate  but  the  driving  function  has 
a  large  amplitude  (  Melvill  Jones  and  Milsum,  1970  ;  Stein,  French 
and  Holden,  1972  )  .  If  there  is  a  carrier  rate  ,  rectification  and 
its  concomitant  distortions  can  be  avoided  by  use  of  small  input 
signals  ,  and  so  in  this  case  rectification  is  a  smooth  static 
nonlinearity  .  However  ,  for  the  case  of  a  zero  carrier  rate 
reduction  of  the  input  amplitude  does  not  lead  to  a  reduction  in  the 
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harmonic  content  of  the  response  ,  and  so  in  this  case  rectification  is 
an  essential  or  ’hard’  nonlinearity  .  Spekreijse  and  Oostings  (  1969  ) 
have  shown  that  the  relative  harmonic  content  may  be  reduced  by  the 
addition  to  tne  driving  function  of  an  auxiliary  signal  which  is  not 
correlated  with  the  driving  function  and  then  averaging  the  response  . 
This  method  of  linearization  by  the  use  of  an  auxiliary  signal  can  be 
used  to  determine  the  order  of  nonlinear  elements  in  a  system  and  has 
been  applied  to  the  analysis  of  visually  evoked  responses  in  man 
(  van  der  Tweel  and  Spekreijse,  1969  )  ,  retinal  ganglion  cell  activity 
(  Spekreijse  ,  1969  )  and  insect  mechanoreceptors  C  Pearson  and  Holden, 
1970  )  . 

Cycle  histograms  of  the  response  of  the  receptor  to  sinusoidal 
driving  functions  added  to  band-limited  white  noise  as  an  auxiliary 
signal  are  compared  in  Figure  4.6  with  the  cycle  histograms  in  response 
to  sinusoidal  driving  functions  .  The  histograms  with  the  auxiliary 
signal  appear  sinusoidal  and  show  neither  the  rectification  nor  the 
phase  locking  apparent  at  the  same  driving  frequency  when  there  is  no 
auxiliary  signal  .  The  auxiliary  signal  ,  which  is  not  correlated 
with  the  driving  function  ,  may  be  considered  to  provide  a  random 
carrier  rate  of  spike  discharge  which  is  modulated  by  the  driving 
function  .  A  Bode  plot  obtained  by  fitting  cycle  histograms  of  the 
response  to  sinusoidal  driving  functions  in  the  presence  of  an 
auxiliary  signal  is  given  in  Figure  4.7  .  The  solid  lines  are  the 
frequency  response  curve  predicted  from  the  step  response  of  Figure 
4.1  .  The  data  points  agree  with  the  predicted  frequency  response 
function  over  a  wider  range  of  frequencies  than  do  the  points  of 
Figure  4.5  .  This  is  because  the  auxiliary  signal  has  broken  the  phase 
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Figure  4.6  .  Cycle  histograms  of  the  response  of  a  tactile  spine  to 

sinusoidal  driving  functions  with  and  without  a  stochastic, 
band-limited  auxiliary  signal  . 
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Figure  4. 
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Bode  plot  of  the  frequency  response  function  obtained 
by  using  sinusoidal  driving  functions  with  a  stochastic  , 
band  ""limited  auxiliary  signal  . 
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locked  pattern  to  gi'/e  a  response  which  when  averaged  is  smoothly 
modulated  . 

Although  use  of  the  auxiliary  signal  has  extended  the  range 
over  which  the  linear  model  of  equation  4.2  applies  it  does  not  give 
any  measure  of  how  accurately  the  model  fits  the  system  . 

Figure  4.8  shows  estimates  of  the  spectra  of  the  tension  driving 
function  (  input  spectrum  )  ,  the  sensory  spike  train  (  output  spectrum  ) , 
and  the  real  (  co-spectrum  )  and  imaginary  (  quad-spectrum  )  parts  of 
the  forward  cross-spectrum  estimated  by  methods  described  in  chapter 
3.  They  are  the  average  of  50  spectral  estimates  ,  each  obtained 
from  2.56  seconds  of  data  sampled  at  5msec.  intervals  .  Thus  each 
spectral  estimate  has  approximately  100  degrees  of  freedom  .  The  input 
spectrum  shows  the  cut-off  characteristics  of  the  filter  used  to  prevent 
aliasing  when  the  tension  record  was  sampled  . 

From  these  spectral  estimates  use  of  equations  3.27  ,  3.28  and 
3.29  provide  an  estimate  of  the  coherence  function  and  its  confidence 
limits  and  equations  3.30  ,  3.32  and  3.33  provide  an  estimate  of  the 
frequency  response  function  and  its  confidence  limits  .  These  estimates 
for  two  receptors  are  shown  in  Figure  4.9  . 

The  extremely  low  level  of  the  coherence  function  indicates  that 
the  system  is  nonlinear  ,  and  so  although  the  frequency  response 
functions  obtained  by  using  different  input  stimuli  (  sine  waves  , 
rectangular  waves  ,  sine  waves  with  auxiliary  functions  and  band- 
limited  white  noise  )  are  all  of  the  same  form  the  frequency  response 
function  is  not  a  good  characterization  of  the  input-output  i.elations 

of  the  receptor  . 
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Figure  4.  8  .  Estimates  of  the  input  ,  output  ,  and  real  and 

imaginary  parts  of  the  cross-spectrum  when  the  tactile 
spine  is  stimulated  by  a  stochastic  ,  band-limited 


driving  function  . 
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Figure  4.  9  .  Bode  plots  of  the  frequency  response  function  ,  and  the 

coherence  function  ,  of  the  tactile  spine  estimated  by 
using  a  stocahstic  ,  band-limited  driving  function  . 
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It  has  been  shown  analytically  that  the  coherence  function  of 
a  perfect  or  linear  voltage  to  rate  convertor  can  be  close  to  one 
at  frequencies  much  less  than  the  carrier  rate  (  Stein,  French  and 
Holden,  1972  )  .  Similar  results  have  been  obtained  using  a  neuronal 
model  of  the  leaky  integrator  type  (  French  and  Stein,  1970  ;  French 

and  Holden,  1971c  ;  Stein,  French  and  Holden,  1972  )  .  Thus  the  low 
value  of  the  coherence  function  of  the  spine  cannot  be  entirely  due 
to  the  threshold  properties  of  the  sensory  neurone  . 

95%  confidence  limits  have  been  drawn  around  the  estimates  of 
the  frequency  response  function  .  These  confidence  limits  depend  on 
the  value  of  the  coherence  function  (  see  section  3.10  )  and  since 
the  coherence  is  low  they  are  very  wide  .  The  narrowing  of  the 
confidence  limits  at  high  frequencies  is  in  part  due  to  the  increase 
in  the  coherence  function  at  higher  frequencies  and  is  also  in  part 
due  to  the  plotting  algorithm  used  .  Since  the  frequency  axis  of 
a  Bode  plot  is  logarithmic  ,  in  order  to  have  the  plotted  points 
approximately  equally  separated  data  points  at  higher  frequencies 
result  from  wider  frequency  segment  averaging  than  points  at  lower 
frequencies  ,  and  so  have  a  higher  number  of  effective  degrees  of 
freedom  . 

From  equation  2.26  the  information  transmission  rate  for  i 
continuous  signals  band-limited  to  50hz.  can  be  estimated  from  the 
coherence  function  estimates  as  6.5  +  1.2  bits/second  for  ten 

preparations  . 
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4.4.  Discussion  . 

The  low  level  of  the  coherence  function  means  that  the  frequency 
response  function  ,  however  it  is  estimated  ,  is  not  a  good  charact¬ 
erization  of  the  input-output  relations  of  the  receptor  .  Bendat 
and  Piersol  (  1966  )  give  three  possible  cases  for  a  coherence 
function  of  a  system  being  less  than  one  : 

a)  .  the  system  is  nonlinear 

b)  .  the  system  is  subject  to  unobserved  inputs  which  are  not 
correlated  with  the  driving  function 

c)  .  the  system  has  an  intrinsic  noise  source  . 

These  three  cases  can  occur  in  any  combination  . 

Since  the  receptor  is  silent  in  the  absence  of  mechanical 
stimulation  of  the  tactile  spine  there  is  no  reason  to  assume  that 
the  low  coherence  of  the  receptor  can  be  attributed  to  unobserved 
inputs  .  Thus  the  low  coherence  can  be  ascribed  to  the  intrinsic 
stochastic  properties  of  the  receptor  and  to  nonlinearites  . 

The  stochastic  properties  of  the  sensory  spike  train  may  not 
be  readily  quantified  as  in  the  absence  of  a  stimulus  there  is  no 
sensory  activity  and  in  the  presence  of  a  maintained  ,  constant 
stimulus  the  spike  train  adapts  and  so  is  non-stationary  . 

A  possible  cause  of  the  low  coherence  is  rectification  due  to 
the  absence  of  a  carrier  rate  .  Figure  4.10  is  a  plot  of  the 
coherence  function  of  a  perfect  rectifier  estimated  by  passing  band- 
limited  white  noise  through  a  Germanium  diode  .  At  no  frequency  is 
the  coherence  function  greater  than  0.75  .  Thus  rectification  alone 
cannot  account  for  the  observed  values  of  the  coherence  function  of 
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Figure  4.10  . 


Estimate  of  the  coherence  function  of  a  Germanium  diode  . 
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the  tactile  spine  . 

One  unusual  aspect  of  rectification  of  a  pulse  train  due  to  the 
absence  of  a  carrier  rate  is  that  the  addition  of  a  DC  or  low  frequency 
signal  incoherent  with  the  input  driving  function  can  increase  the 

coherence  at  all  frequencies  .  Figure  4.11  shows  the  coherence 

\ 

function  of  a  model  neurone  (  French  and  Stein,  1970  )  estimated 
with  no  bias  on  the  integrator  and  so  a  zero  carrier  rate  and  with 
a  bias  on  the  integrator  producung  a  carrier  rate  of  100  impulses/sec.  . 
The  coherence  function  was  estimated  by  subjecting  the  model  neurone 
to  band-limited  white  noise  which  in  both  cases  had  an  input  signal 
root  mean  squared  power  of  0.25  volts*  ‘T.he  time  constant  of  the 
leaky  integrator  was  10  msecs.  .  The  addition  of  a  DC  input  to  the 
integrator  has  increased  the  coherence  at  all  frequencies  ,  and  so 
has  increased  the  information  transmission  rate  estimated  by  equation’ 
2.26  by  30%  .  In  a  linear  system  the  addition  of  any  signal  uncorrelated 
with  the  input  would  cause  a  decrease  in  the  coherence  of  the  system 
at  some  frequencies  and  would  not  increase  the  coherence  at  any 
frequencies  . 

Crowe  (  1967  )  has  pointed  out  the  similarity  of  the  ramp  response 
of  the  tactile  spine  to  the  ramp  response  of  mammalian  muscle  spindle 
secondary  afferents  .  The  coherence  function  of  the  amphibian  muscle 
spindle  afferent  has  been  estimated  by  a  different  technique  (  Koles, 
1970  )  and  although  the  values  of  the  coherence  function  of  the 
spindle  are  greater  than  those  of  the  tactile  spine  the  coherence 
function  has  the  same  general  shape  ,  that  is  ,  increasing  at  higher 
frequencies  .  This  contrasts  with  the  coherence  function  computed 
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Figure  4.  11  .  Estimates  of  the  coherence  function  of  a  model  neurone 

with  a  zero  carrier  rate  and  a  carrier  rate  of  100 
impulses/second  caused  by  an  auxiliary  DC  input  . 
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(  Stein,  French  and  Holden, 1972  )  or  estimated  (  French  and  Holden, 
1971c  )  for  simple  neuronal  models  ,  for  which  the  coherence  function 
decreases  with  frequency  .  This  increase  in  the  coherence  function 
with  frequency  may  be  due  to  the  processes  underlying  adaptation  . 

The  effect  of  adaptation  (  produced  by  increasing  the  absolute 
refractory  period  )  on  the  coherence  function  of  the  model  neurone 
is  to  reduce  the  coherence  at  low  frequencies  while  leaving  the 
coherence  at  high  frequencies  uneffected  .  This  supports  the  idea 
that  adaptation  gives  a  receptor  a  'high-pass’  filter  characteristic  . 

Even  though  the  frequency  response  function  is  not  a  good 
characterization  of  the  input-output  relations  of  the  tactile  spine 
the  estimates  obtained  by  using  spectral  analysis  may  be  fitted 

i 

by  the  frequency  response  function  of  equation  4.2  .  An  expression 
of  this  form  ,  derived  from  the  power  function  decay  in  response  to 
a  step  input  function  ,  has  been  used  to  model  the  input-output 
relations  of  a  wide  range  of  receptors  -  the  tactile  spine  of  the 
cockroach  (  Chapman  and  Smith,  1963  ;  Pearson  and  Holden,  1970  )  , 
the  slowly  adapting  stretch  receptor  of  the  crayfish  (  Brown  and 
Stein,  1966  )  ,  carotid  body  chemoreceptors  (  Landgren,  1952  )  and 
abdominal  stretch  receptors  in  butterflies  (  Weevers,  1966  )  .  Various 
physical  analogs  have  been  proposed  which  give  a  power  function 
decay  in  response  to  a  step  input  ,  and  a  simple  abstract  model  is 
a  series  of  exponential  decay  processes  which  have  different  time 
courses  . 

Chapman  and  Smith  (  1963  )  describe  the  power  function  decay 
as  a  fractional  differentiation  of  the  driving  signal  .  Brown  and 
Stein  (  1966  )  point  out  that  the  power  function  provides  a  means  of 
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encoding  information  about  the  driving  function  and  its  first 
derivative  .  However  ,  because  of  the  low  coherence  of  the  receptor  , 
an  alternative  view  would  be  to  consider  it  as  a  detector  of  any 
sharp  or  sudden  sudden  stimulus  -  although  the  waveform  of  the 
stimulus  cannot  be  reliably  recovered  from  the  spike  train  the  fact 
that  a  transient  had  occured  could  be  detected  . 


/ 


1  ••  u  '  1  •  ■  ^  g  l  ^  '  " 


71 


CHAPTER  5  . 

THE  FEMORAL  CHORDOTONAL  ORGAN  OF  THE  LOCU ST  . 

5.1  .  Introduction. 

Insect  mechanoreceptors  and  proprioceptors  have  been  divided  into 
five  morphological  groups  :  hair  plates  ,  campaniform  sensilla  ,_stretch 
receptors  ,  chordotonal  organs  and  statocys t-like  organs  (  Dethier, 

1963  )  .  The  chordotonal  organs  are  morphologically  the  most  complex 
kind  of  mechanoreceptor  found  in  insects  .  A  chordotonal  organ 
consists  of  a  number  of  sensilla  attached  to  the  body  wall  by  a 
ligament  at  the  axonal  end  and  by  accessory  cells  at  the  dendritic 
end  .  Each  sensillum  consists  of  a  bipolar  sensory  neurone  and  at 
least  two  accessory  cells  .  The  dendrite  of  the  sensory  neurone  is 
enveloped  by  the  two  or  more  accessory  cells  ,  and  its  termination 
is  enclosed  by  a  scolopoid  sheath  .  In  spite  of  extensive  electro- 
physiological  and  histological  investigations  of  Crustacean  chordo¬ 
tonal  organs  (  Wiersma  and  Boettiger,  1959  ;  Whitear,  1962  ; 

Mendelson,  1963  ;  Bush,  1965a  ,  1965b  )  the  significance  of  the 
complex  structure  of  the  chordotonal  organ  is  unknown  . 

Simple  ,  connective  chordotonal  organs  (  Howse,  1968  )  in  insects 
have  not  been  investigated  in  any  great  detail  .  Becht  (  1958  ) 
has  described  the  effects  of  insecticides  on  cockroach  chordotonal 
organs  ,  and  Hubbard  (  1959  )  has  published  a  preliminary  study 
of  the  femoral  chordotonal  organ  of  the  locust  . 
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As  part  of  a  series  of  experiments  investigating  the  role  of 
proprioceptive  feedback  during  postural  and  locomotor  activity  in 
locusts  Runion  and  Usherwood  (  1966  )  and  Usherwood,  Runion  and 
Campbell  (  1968  )  have  described  the  histology  of  the  femoral 
chordotonal  organ  of  the  locust  and  have  qualitatively  investigated 
its  responses  and  reflex  actions  . 

The  femoral  chordotonal  organ  is  attached  by  ligaments  to  the 
posterior  wall  of  the  femur  and  by  an  apodeme  to  a  protrusion  on  the 
tibia  near  the  f emoro-tibial  articulation  ,  It  is  also  attached  by 
ligaments  to  the  apodemes  of  the  extensor  and  flexor  tibiae  muscles  . 
Three  groups  of  sensory  cells  have  been  described  ,  comprising  of  a 
total  of  24  bipolar  sensory  neurones  in  the  chordotonal  organ  .  The 
sensory  axons  from  the  chordotonal  organ  are  contained  in  metathoracic 
nerve  N<-b^  [  terminology  of  Campbell,  1961  ]  . 

Recordings  from  the  chordotonal  organ  nerve  in  response  to 
sinusoidal  or  ramp  changes  in  the  f emoro-tibial  angle  showed  that 
there  were  both  tonic  and  phasic  units  .  It  was  not  possible  to 
activate  single  units  ,  but  the  summed  tonic  discharge  was  a  linear 
function  of  the  f emoro-tibial  angle  .  The  summed  tonic  discharge 
adapted  to  within  90%  of  a  steady  level  within  one  minute  after  a 
step  change  in  f emoro-tibial  angle  . 

Phasic  units  responded  to  sinusoidal  movements  of  frequencies 
from  0.1  to  170  hz.  ,  some  units  responding  to  extension  and  others  to 
flexion  .  Chordotonal  sensory  activity  during  extension  of  the 
tibia  caused  reflex  activation  of  the  flexor  tibiae  muscle  ,  and  during 
flexion  caused  reflex  activation  of  the  slow  axon  in  nerve  N^b  to 
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the  extensor  tibiae  muscle  .  Generally  the  response  to  extension 
was  greater  than  the  response  to  flexion  .  With  one  hz.  sinusoidal 
stimulation  the  reflex  response  to  the  extensor  muscle  adapted  . 

5.2  .  Methods  . 

In  all  experiments  adult  male  locusts  (  Scliistocerca  gregaria  )  , 
which  had  been  commercially  obtained  ,  were  used  .  The  metathoracic 
leg  was  cut  through  at  the  coxa  ,  and  the  cut  surface  covered  with 
Vaseline  to  prevent  drying  .  The  leg  was  mounted  on  a  cork  platform 
so  that  the  tibia  was  held  in  the  vertical  plane  .  The  leg  was  held 
rigidly  by  dental  wax  ,  plasticine  or  glue  (  Eastman  9-10  adhesive  )  . 
The  tibia  was  cut  5  mm.  from  the  f emoro-tibial  joint  and  tissue 
inside  the  tibial  stump  macerated  by  a  blunt  probe  .  The  femoro- 
tibial  joint  was  dissected  so  that  the  tibia  was  attached  to  the 
femur  only  by  the  apodemes  of  the  chordotonal  organ  and  the  flexor 
and  extensor  tibiae  muscles  .  The  tibial  stump  was  held  horizontally 
in  a  pair  of  screw  forceps  attached  to  the  armature  of  a  vibration 
generator  (  Pye-Ling  model  V47  )  .  A  small  section  of  the  ventral 
cuticle  of  the  femur  two  cms .  from  the  f emoro-tibial  joint  was 
removed  and  the  underlying  flexor  tibiae  muscle  dissected  away  ,  and 
picked  up  on  bipolar  silver  electrodes  .  The  f emoro-tibial 
joint  and  the  recording  electrodes  were  covered  with  Vaseline  to 

i 

prevent  drying  .  This  minimal  dissection  avoids  exposing  the 
chordotonal  organ  ,  as  Usherwood,  Runion  and  Campbell  (  1968  )  have 
observed  that  Vaseline  or  locust  saline  (  Hoyle,  1953  )  change  the 
responsiveness  of  the  chordotonal  organ  . 
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Movement  of  the  f emoro-tibial  joint  is  limited  to  an  angle  of 
170  degrees  ,  which  corresponds  to  a  1.9  mm.  movement  of  the 
apodeme  of  the  chordotonal  organ  .  All  movements  of  the  tibia  were 
restricted  within  _+  0.5  mm.  of  the  extension  equivalent  to  90 
degrees  flexion  .  The  vibration  generator  was  driven  by  the  output 
from  a  length  servo  control  unit  via  a  power  amplifier  ,  The 
feedback  transducer  to  the  servo  was  a  Hewlett-Packard  DC  input, 

DC  output  length  transducer  (  model  24DCDT  050  )  which  had  a  linear 
displacement  range  of  +  0.050  inches  .  The  command  signal  to  the 
servo  could  be  a  DC  level  added  to 

i) .  the  output  from  a  Hewlett-Packard  3300A  function  generator 
ii) .  the  output  from  a  band-limited  noise  generator  . 

The  closed  loop  frequency  response  function  of  the  length  servo 
control  system  was  flat  up  to  100  hz .  . 

The  amplified  action  potentials  ,  the  output  from  the  length 
transducer  and  reference  markers  were  recorded  on  a  Thermionic 
Products  T3000  FM  tape  recorder  for  subsequent  processing  on  the  LAB-8 
computer  .  Data  anlysis  was  by  methods  described  in  chapters  3  and  4  . 

5 .3  .  Results  . 

An  interspike  interval  histogram  of  the  discharge  of  a  tonic 
unit  from  the  chordotonal  organ  is  shown  in  Figure  5.1a  .  Since  an 
increase  in  the  length  of  the  chordotonal  organ  ,  corresponding  to 
an  increase  in  the  f emoro-tibial  angle  ,  causes  recruitment  of  other 
tonic  units  records  in  which  a  single  tonic  unit  could  be  reliably 
extracted  by  pulse  height  analysis  could  only  be  obtained  with  low 
discherge  rates  .  The  unit  of  Figure  5.1a  had  a  mean  rate  of  10.3 
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Figure  5.1 


A 


100  MSECS 


Interspike  interval  histograms  of  two  tonic  units  from 
the  femoral  chordotonal  organ  .  Inserts  show  computed 


Gamma  distributions  . 
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impul ses/ second  and  a  coefficient  of  variation  of  the  interspike 
interval  of  0.63  ..  Since  the  interspike  interval  histogram  is 
unimodal  and  positively  skewed  it  seems  appropriate  to  fit  it  by  a 
Gamma  distribution  (  Cox  and  Lewis,  1966  )  .  The  Gamma  distribution 
is  defined  by  : 


f  (x) 


r  k  ,k  k-1  -kx/y 

[  u  ]  -* - § - 

r  (k) 


5.1 


where  y-  E  {  x}  and  T  (k)  is  the  Gamma  function  of  order  k 
defined  by  equation  4.3  .  The  variance  of  this  distribution  is  y  /k 
and  the  coefficient  of  variation  is  1//  k  .  When  k  is  one  the 
Gamma  distribution  reduces  to  the  exponential  distribution  ,  and  when 
k  tends  to  infinity  the  Gamma  distribution  tends  to  the  Normal 
distribution  .  Thus  the  Gamma  distribution  can  encompass  a  wide 
range  of  distributions  ,  however  ,  it  does  not  encompass  any  negatively 
skewed  ,  flat  ,  sharply  peaked  or  multimodal  distributions  . 

k  ,  the  order  of  the  Gamma  distribution  may  be  estimated  from  : 

1 

k  'v  - - -  5.2 

t  -  i  2 

V 

where  a  is  the  standard  deviation  .  Use  of  equation  5.2  with 


estimates  of  the  mean  and  standard  deviation  of  the  histogram  of 
Figure  5.1a  give  the  estimate  of  k  as  2.7  .  The  insert  of  Figure 
5.1a  shows  the  computed  Gamma  distribution  of  order  2.7. 

This  method  of  estimating  k  does  not  produce  the  least  mean 

A  A 

square  error  between  k  and  k  ,  especially  as  k  <  10  (  Cox  and 
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Lewis,  1966)  and  oO  is  not  an  efficient  method  of  estimating  k  (  Bendat 
and  Piersol ,  1966  )  .  Cox  and  Lewis  (  1966  )  have  shown  that  a 
maximum  liklihood  estimator  of  k  is  given  by  the  solution  of 

A 

log  (k)  -  Y  (k)  =  log  t  -  "  log  t. 
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and  T  (k)  is  the  digamma  function  defined  by 
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Some  tonic  units  fired  regularly  :  Figure  5.1b  shows  the 
interspike  interval  histogram  of  such  a  unit  ,  with  a  mean  rate'  of 
20.1  impulses/second  and  a  coefficient  of  variation  of  0.24  . 

From  equation  5.2  this  would  give  an  estimate  of  the  order  of  the 
distribution  k  =  17.3  .  The  computed  Gamma  distribution  of  order 

17.3  is  shown  in  the  insert  of  Figure  5.1b  .  The  power  spectrum 
estimate  of  this  spike  train  obtained  from  50  non-overlapping  samples 
each  of  5.12  seconds  is  given  in  Figure  5.2  .  This  spectrum  is 
typical  of  a  high  order  Gamma  distribution  ,  showing 
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the  low  and  high  frequency  asymptotes  and  the  overshoot  at  the  mean 
rate  . 

Cycle  histograms  of  the  response  of  a  single  unit  to 
sinusoidal  driving  functions  are  shown  in  Figure  5.3  .  In  contrast 
to  cycle  histograms  obtained  from  the  tactile  spine  ,  the  histograms 
obtained  from  the  chordotonal  organ  do  not  show  rectification  at  low 
driving  frequencies  ,  This  is  because  at  a  constant  length  the 
receptor  has  a  variable  background  discharge  or  carrier  rate  which 
may  be  smoothly  modulated  by  low  frequency  sinusoidal  driving 
functions  .  In  the  case  of  the  tactile  spine  of  the  cockroach 
rectification  was  proposed  as  a  partial  cause  of  the  low  coherence 
function  observed  .  Thus  ,  if  other  factors  were  equal  ,  one  might 
expect  that  the  coherence  function  of  the  chordotonal  organ  would 
be  greater  than  the  coherence  function  of  the  tactile  spine  . 

Phase  locking  in  response  to  sinusoidal  driving  functions  occurs 
at  a  frequency  higher  than  that  for  the  tactile  spine  .  At  50  hz . 
there  is  one-to-one  entrainment  » 

Figure  5.4  shows  estimates  of  the  spectra  of  the  input  driving 
function  (  input  spectrum  )  ,  the  output  spike  train  of  a  single  unit 
(  output  spectrum  )  ,  the  real  part  (  co-spectrum  )  and  imaginary 
part  (  quad-spectrum  )  of  the  cross-spectrum  when  the  chordotonal 
organ  was  excited  by  a  stochastic  ,  band-limited  driving  function  . 
These  estimates  are  the  average  of  50  sample  estimates  each  from  5.12 
seconds  of  data  .  All  the  estimates  have  been  normalized  to  the 
eleven  bit  accuracy  of  SAPA  (  French  and  Holden,  1971b  )  and  so  no 
amplitude  scale  is  given  .  Since  there  has  been  no  frequency  range 
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Figure  5.3  .  Cycle  histograms  of  the  response  of  a  phasic  unit  from 
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Figure  5.4  .  Estimates  of  the  input  ,  output  and  cross-spectra  when 

the  chordo tonal  organ  was  excited  by  a  stochastic  , 
band-limited  driving  function  .  The  ordinates  are  a 
logarithmic  plot  of  the  SAPB  file  contents  . 
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averaging  the  estimates  have  approximately  100  degrees  of  freedom  . 

Estimates  of  the  frequency  response  function  and  its  confidence 
limits  obtained  by  using  equations  3.30,  3.32  and  3.33  and  of  the 
coherence  function  obtained  by  using  equations  3.27,  3.28  and  3,29' 
are  shown  in  figure  5.5  .  Because  of  the  low  value  of  the  coherence 
function  the  95%  confidence  limits  of  the  sensitivity  curve  are  wide  , 
but  for  frequencies  less  than  10  hz.  the  sensitivity  is  approximately 
constant  .  Above  10  hz.  the  sensitivity  increases  .  This  simple 
frequency  response  curve  ,  illustrated  for  four  different  units  in 
Figure  5.6  ,  is  similar  to  that  of  a  simple  high  pass  filter  whose 
transfer  function  is  given  by 

s 

H(s)  =  (  sT  +  1  )  5.6 


A  Bode  plot  of  this  transfer  function  is  shown  in  Figure  5.7  .  The 
corner  frequency  of  the  gain  curve  occurs  at  1/T  hz .  ,  above  which 
the  gain  increases  at  +  6  dB  /  octave  .  The  slope  of  the  frequency 
response  curves  of  Figure  5.5  above  10  hz.  is  about  100  impulses/ 
second,  mm  /  decade  change  in  frequency  .  Since  below  10  hz.  the 
sensitivity  is  approximately  constant  at  10  impulses/sec.  mm.  if 
a  sensitivity  of  10  impulses/sec.  mm.  is  taken  as  equivalent  to  a 
gain  of  one  the  slope  of  the  frequency  response  curve  above  10  hz. 
is  equivalent  to  +20  dB  /  decade  ,  or  +6  dB  /  octave  .  This  arguement 
is  not  strictly  valid  since  the  input  and  output  do  not  have  the 
same  dimensions  ,  and  so  a  gain  may  not  be  defined  ;  however  ,it 
does  support  the  simple  frequency  response  function  of  equation  5.6  . 
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Figure  5.5  .  Estimates  of  the  frequency  response  function  of  two 

chordotonal  organ  units  .  The  sensitivity  curves  are 
in  impulses/sec .  mm.  . 
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Figure  5.6  . 


Sensitivity  of  four  chordotonal  organ  units  . 
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Bode  plot  of  trans 


Figure  5.7  . 
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Even  though  the  .95%  confidence  limits  are  wide  the  pha,se  advance 

of  the  frequency  response  functions  of  Figure  5.5  do  support  the  simple 

frequency  response  function  of  equation  5.6  in  that  the  phase  advances 

monotonically  from  zero  to  +  90  degrees  .  However  because  of  the  wide 

* 

confidence  limits  local  maxima  suggestive  of  other  terms  in  the 
frequency  response  function  would  not  be  detected 

Figure  5.8  shows  estimates  of  the  coherence  function  of  four 
different  units  .  The  coherence  function  is  less  than  0.5  at  all 
frequencies  ,  and  extremely  low  at  frequencies  less  than  10  hz .  These 
coherence  function  curves  have  a  similar  shape  to  those  for  the 
cockroach  tactile  spine  (  Chapter  4  )  and  the  amphibian  muscle  spindle 
(  Koles,  1970  )  .  The  value  of  the  coherence  function  is  higher  than 
that  for  the  cockroach  tactile  spine  ,  as  was  predicted  from  the 
absence  of  rectification  in  the  cycle  histograms  in  response  to 
sinusoidal  driving  functions  . 

The  arrows  on  the  coherence  function  curves  of  Figure  5.8  mark 
the  mean  firing  rate  of  the  unit  .  In  experiments  using  an  electronic 
neural  analog  (  French  and  Stein,  1970  )  the  coherence  function  drops 
sharply  before  the  mean  firing  rate  or  carrier  rate  and  is  low  above 
the  carrier  rate  (  French  and  Holden,  1971c  ;  Stein,  French  and 
Holden,  1972  )  .  This  is  not  observed  in  the  coherence  function 
estimates  for  the  chordo tonal  organ  .  This  can  be  explained  as  in 
the  case  of  the  neural  anlog  the  driving  function  modulated  the  carrier 
rate  ,  which  was  the  mean  firing  rate  .  However  ,  in  the  case  of  the 
chordotonal  organ  ,  the  resting  rate  is  variable  ,  and  the  mean  rate 
is  caused  by  the  phasic  response  to  the  driving  function  ,  a.nd  so  is 
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Figure  5.8  .  Coherence  function  of  four  chordotonal  organ  units  . 
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not  entirely  analogous  to  a  carrier  rate  . 

The  coherence  function  of  the  electronic  neural  analog  is  a 
function  of  the  input  RMS  power  (  French  and  Holden,  1971c  ;  Stein, 
French  and  Holden,  1972  )  .  Figure  5.9  shows  coherence  function 
estimates  for  a  single  chordotonal  organ  unit  at  three  different  RMS 
power  values  of  input  length  driving  function  .  The  estimates  have 
more  than  a  hundred  degrees  of  freedom  .  The  coherence  increases  with 
the  input  amplitude  ,  especially  at  high  frequencies  .  This  contrasts 
with  the  results  from  the  electronic  neural  analog  ,  in  which  the 
increase  in  coherence  with  input  amplitude  was  predominantly  at 
frequencies  less  than  the  carrier  rate  .  This  increase  in  coherence 
at  high  frequencies  supports  the  idea  that  the  receptor  acts  as  a 
high  pass  filter  . 

The  information  transmission  rate  of  chordotonal  organ  single 
units  calculated  by  equation  2.26  ranged  from  8.75  to  27.7  bits/second  . 
Since  the  coherence  function  increases  with  input  amplitude  the 
information  transmission  rate  increases  with  the  amplitude  of  the 
driving  function  . 

5.4  .  Discussion  . 

The  frequency  response  functions  obtained  by  spectral  analysis 
indicate  that  although  the  coherence  function  is  low  ,  and  so  a  linear 
model  is  not  an  accurate  characterization  of  the  input-output  relations 
of  the  receptor  ,  the  best- fitting  linear  model  is  very  similar  to 
a  simple  lead  filter  whose  transfer  function  is  given  by  equation 

5.6  and  shown  in  Figure  5.7  .  A  possible  cause  for  this  extremely 
simple  transfer  function  might  be  the  visco-elastic  properties  of  the 
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Figure  5.9  .  Estimates  of  the  coherence  function  of  a  chordotonal 

organ  unit  at  three  different  input  amplitudes  . 
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apodeme  of  the  chordotonal  organ  , 

The  low  coherence  of  chordotonal  organ  single  units  at  frequencies 
below  10  hz.  means  that  the  waveform  of  slow  changes  of  the  femoro— 
tibial  angle  may  not  be  accurately  recovered  from  the  spike  train  of 
a  single  unit  by  any  linear  operation  .  Changes  in  the  femoro- 
tibial  angle  at  frequencies  from  10  hz.  to  50  hz.  may  only  be 
partially  recovered  by  a  linear  operation  on  a  single  spike  train  . 

The  dependence  of  the  coherence  function  on  the  amplitude  of  the 
driving  function  means  that  more  information  may  be  recovered  by 
a  linear  operation  on  the  spike  train  about  large  amplitude 
movements  than  about  small  amplitude  movements  . 

Usherwood ,Runion  and  Campbell  (  1968  )  have  reported  that 
destruction  of  the  femoral  chordotonal  organs  produce  significant 
changes  in  postural  and  locomotor  behaviour  .  This  has  been  used 
to  support  the  concept  that  proprioceptive  feedback  is  very 
important  in  the  generation  of  normal  postural  and  walking 
behaviour  .  This  concept  implies  that  the  proprioceptors  can 
provide  the  central  nervous  system  with  detailed  information  about 
position  and  movement  ,  and  that  this  information  can  be  ’decoded* 
or  utilized  by  the  central  nervous  system  in  the  normal  control 
of  movement  .  During  walking  and  postural  activity  movements  of  the 
f emoro- t ib ial  angle  are  slow  ,  and  occur  at  frequencies  where  the 
coherence  of  single  chordotonal  organ  units  is  very  low  .  This  low 
coherence  means  that  detailed  proprioceptive  information  cannot  be 
recovered  from  the  spike  train  of  a  chordotonal  organ  single  unit 
by  any  linear  operation  .  Mountcastle  (  1967  )  has  reviewed 


V 


"  •  ‘  ■ 


91 


evidence  that  supports  the  hypothesis  that  the  processing  of  sensory 
information  by  the  nervous  system  is  by  approximately  linear  operations  . 
If  the  decoding'  operations  performed  on  the  spike  train  of  a  single 
chordotonal  organ  unit  are  approximately  linear  then  the  single  unit 
cannot  provide  the  central  nervous  system  with  detailed  proprioceptive 
information  .  This  could  be  used  to  argue  against  the  importance  of 
proprioceptive  feedback  from  the  chordotonal  organ  during  the  normal 
control  of  postural  and  locomotor  activity  . 

However  ,  each  femoral  chordotonal  organ  contains  about  24  bipolar 
sensory  neurones  .  If  the  central  nervous  system  averaged  the  response 
of  all  those  chordotonal  organ  units  which  had  the  same  properties 
then  the  chordotonal  organ  might  provide  the  central  nervous  system 
with  detailed  proprioceptive  information  which  could  be  recovered  by 
a  linear  operation  . 

This  would  be  supported  if  the  coherence  between  the  length  of 
the  chordotonal  organ  and  the  reflex  motor  discharge  was  greater 
than  the  coherence  between  the  length  of  the  chordotonal  organ  and  a 
single  unit  .  However  ,  preliminary  experiments  with  headless  locusts 
give  a  very  low  coherence  between  length  of  the  chordotonal  organ  and 
the  reflex  discharge  in  the  slow  axon  of  N^b  to  the  extensor  tibiae 
muscle  .  This  low  coherence  may  be  due  to  other  inputs  to  the  motor 
neurone  which  might  be  removed  by  isolating  the  ganglion  from  all 
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CHAPTER  6  . 


DISCUSSION  . 

All  physical  systems  are  nonlinear  ,  and  so  if  linear  systems 
theory  is  to  be  used  to  formulate  an  abstract  model  to  describe  the 
input-output  relations  of  a  physical  system  some  measure  of  the 
accuracy  of  the  model  is  essential  .  When  the  physical  system  is  a 
sensory  neurone  which  responds  to  a  continuous  driving  function  by 
a  train  of  action  potentials  there  is  considerable  non— harmonic 
distortion  ,  and  so  harmonic  distortion,  however  measured  ,  is  not 
a  good  measure  of  the  linearity  of  the  system  .  The  coherence  function  , 
defined  by  equation  2.22  ,  provides  a  measure  of  the  linearity  of  the 
relation  between  the  input  and  output  which  is  sensitive  to  all  kinds 
of  nonlinearities  and  to  intrinsic  noise  .  The  coherence  function 
can  only  be  estimated  by  subjecting  the  system  to  a  broad  band  driving 
function  ,  and  using  consistent  ,  alias  and  bias  free  spectral  analysis  . 
If  that  part  of  the  output  which  is  linearly  related  to  the  input 
driving  function  is  considered  to  be  the  signal  ,  and  that  part  of 
the  output  which  is  not  linearly  related  to  the  input  is  considered  to 
be  noise  ,  then  the  coherence  function  is  related  to  the  signal-to- 
noise  ratio  as  shown  in  equation  2.24  .  If  the  input  driving 
function  is  a  sinusoid  ,  which  has  a  line  spectrum  ,  then  the  signal- 
to-noise  ratio  will  be  infinite  and  so  the  coherence  function  will  be 
equal  to  one  .  Thus  inputs  which  have  line  or  very  narrow  band  spectra  , 
which  are  simple  deterministic  functions-  of  time  ,  cannot  Be  used  to 
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estimate  the  coherence  function  , 

There  are  a  variety  of  methods  which  may  Be  used  to  estimate  the 
spectral  density  matrix  of  section)  2.3  .  The  rediscovery  of  the 
Fast  Fourier  Transform  algorithm  (  Cooley  and  Tukey,  1965  )  has  made 
the  direct  method  of  spectral  estimation  computationally  economical  , 
and  the  development  of  alias-free  sampling  of  stochastic  point  processes 
with  a  regular  sampling  interval  (  French  and  Holden,  1971a  ,  1972  ) 
has  permitted  the  use  of  the  Fast  Fourier  Transform  algorithm  in  the 
estimation  of  the  spectrum  of  a  spike  train  [  when  the  spike  train  is 

’i 

treated  as  a  realization  of  a  stochastic  point  process  or  as  a  series 
of  unit  impulses  ]  and  of  the  cross-spectrum  between  two  spike  trains 
or  between  a  spike  train  and  a  continuous  signal  . 

The  direct  method  of  spectral  estimation  has  been  implemented 
on  a  small  computer  ,  and  the  consistency  of  the  estimates  has  been 
achieved  by  averaging  spectral  estimates  from  non-overlapping  samples  . 

A  better  method  of  consistent  spectral  estimation  would  be  to 
average  spectral  estimates  from  samples  with  a  50%  overlap  (  Welch, 

.1967  )  .  However  ,  implementation  of  this  technique  on  a  small 
computer  would  necessitate  a  large  ,  bulk  storage  device  such  as  a 
magnetic  tape  unit  in  order  to  store  the  entire  record  and  then 
sequentially  access  overlapping  segments  .  If  such  a  storage  device 
were  available  it  would  permit  the  computation  of  spectral  estimates 
with  the  number  of  data  points  greater  than  the  number  of  core  memory 
locations  (  Brenner,  1969  )  . 

Consistent  spectral  estimates  provide  a  means  of  estimating  the 
frequency  response  function  as  well  as  the  coherence  function  ,  and 
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from  the  coherence  function  and  the  number  of  degrees  of  freedom  of 
the  spectral  estimates  the  confidence  limits  of  the  frequency  response 
function  can  be  computed  .  In  both  the  tactile  spine  of  the  cockroach 
and  the  chordotonal  organ  of  the  locust  the  sensitivity  of  the 
frequency  response  function  increased  at  higher  frequencies  ,  as 
did  the  coherence  .  Thus  the  estimate  of  the  frequency  response  function 
was  less  stable  at  low  frequencies  than  at  high  frequencies  .  If  the 
input  driving  function  had  been  shaped  so  that  there  was  more  power 
present  at  low  frequencies  than  at  high  frequencies  the  stability  of 
the  frequency  response  function  would  have  been  more  even  . 

From  the  coherence  function  estimate  the  information  transmission 
rate  of  the  receptor  could  be  estimated  by  use  of  equation  2.26  . 

This  is  not  an  estimate  of  the  information  capacity  of  the  receptor 
as  the  information  capacity  is  defined  as  the  maximum  possible  transmission 
rate  .  This  would  only  be  achieved  if  the  input  driving  function  had 
a  spectrum  so  that  most  of  the  input  power  was  at  frequecies  where 
the  coherence  was  highest  .  Previous  estimates  of  the  information 
transmission  rate  of  receptors  have  been  limited  to  quantifying  the 
ability  of  the  receptor  to  transmit  information  about  a  steady  level 
(  Werner  and  Mountcastle,  1965  ;  Darian- Smith  et  al,  1968  ;  Stein, 

1968b  ;  Matthews  and  Stein,  1969b  )  .  Thus  the  estimates  obtained  by 
use  of  equation  2.26  for  the  tactile  spine  and  the  chordotonal  organ 
are  higher  than  those  reported  for  other  receptors  . 

Since  the  nervous  system  is  made  up  of  a  large  number  of  neurones 
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it  might  Be  argued  that  the  detailed  analysis  of  the  input-output 
relations  of  a  single  receptor  does  not  help  very  much  In  understanding 
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the  functioning  of  the  nervous  system  .  The  spectral  density  matrix 
of  section  2.3  may  be  generalized  for  a  system  with  an  arbitrary 
number  of  inputs  and  outputs  .  Thus  the  methods  described  in 
Cnapters  2  and  3  and  applied  in  Chapters  4  and  5  can  be  generalized  to 
systems  with  several  inputs  and  several  outputs  .  The  complexity  of 
the  system  to  be  analyzed  by  these  methods  is  only  restricted  by  the 
number  of  units  which  may  be  recorded  and  the  size  of  the  data  processing 
system  .  Thus  when  it  is  possible  to  record  from  a  large  number 
of  neurones  simultaneously  the  methods  of  Chapters  2  and  3  could 
provide  a  means  of  dealing  with  the  large  amount  of  data  . 

The  driving  function  to  the  receptors  studied  was  a  function  of 
time  only  .  However  ,  many  sensory  systems  respond  to  driving  functions 
which  are  functions  of  spatial  variables  as  well  as  of  time  .  Spectral 
analysis  is  readily  generalized  to  include  n-dimensional  functions 
and  so  the  methods  developed  in  this  thesis  can  provide  the 
background  to  an  approach  to  more  complex  aspects  of  sensory 


information  transmission  . 
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